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I.  INTRODUCTION 


Digital  signal  processing  is  a  field  which  is  rapidly  expanding  due  to  advances  in 
modern  technology.  Essential  to  this  field  are  digital  filters.  Modeling  these  filters  con¬ 
stitutes  much  of  the  effort  involved  in  digital  signal  processing.  The  filters  provide  a 
transfer  function  which  describes  the  relationship  between  filter  input  and  output. 
Autoregressive  (AR),  moving  average  (MA)  and  combination  autoregressive  moving 
average  (ARMA)  models  are  widely  used  to  represent  the  transfer  function  of  a  digital 
filter.  A  filter  transfer  function  is  commonly  described  in  direct  form.  This  form  is  a  ratio 
of  two  polynomials,  usually  of  the  form, 


H{z)  = 


B(z) 

A(z) 


b^  4-  b\  z  +•••-+ '  b(z 

1  +  ax  z~]  A —  +  as  z~s 


(1.1) 


The  above  equation  describes  an  ARMA  model  of  order  (s,t)  where  s  is  the  order  of  the 
denominator  and  t  the  order  of  the  numerator.  The  a parameters  form  the 
autoregressive  portion  of  the  ARMA  model.  The  b,  parameters  form  the  moving  average 
portion  of  the  ARMA  model.  If  all  the  autoregressive  parameters  are  zero,  then  the  filter 
transfer  function  H(z)  is  strictly  a  moving  average  process  of  order  t.  If  all  the  moving 
average  parameters  are  zero  except  for  b0  equal  to  one,  then  the  filter  transfer  function 
is  strictly  an  autoregressive  process  of  order  s. 

A.  OBJECTIVES  OF  THE  THESIS 

Fundamental  to  the  design  of  digital  filters  is  estimation  of  AR,  MA  or  ARMA 
parameters.  Accurate  and  efficient  parameter  estimation  has  been  the  subject  in  much 
of  the  related  digital  signal  processing  literature  [Refs.  1  ,2,3].  The  first  objective  of  this 
thesis  is  to  confirm  the  proposed  ARMA  parameter  estimation  algorithm  of  [Ref.  A:  pp. 
619-621],  which  leads  to  the  design  of  a  new  ARMA  digital  lattice  filter.  The  proposed 
algorithm  is  a  generalization  of  the  Mullis-Roberts  criterion  for  parameter  estimation 
known  as  the  modified  least  squares  problem  [Ref.  5:  pp.  227-228].  The  algorithm  uses 
two  recursive  formula  to  estimate  the  parameters.  One  is  an  A R  recursive  formula  which 
estimates  ARMA  parametes  as  the  AR  order  is  increased  by  one.  The  other  is  an  MA 
recursive  formula  which  estimates  ARMA  parameters  as  the  MA  order  is  increased  by 
one.  This  algorithm  is  unique  in  that  it  allows  for  the  design  of  an  ARMA  model  with 
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arbitrary  AR  and  MA  orders  with  no  dependency  of  an  AR  model  on  an  MA  model  or 
vice  versa.  The  ARMA  digital  filter  designed  from  the  proposed  ARMA  parameter  es¬ 
timation  algorithm  is  in  the  form  of  a  lattice  structure.  Lattice  realizations  of  filters  are 
widely  used  and  provide  excellent  analysis  of  prediction  errors  [Refs.  6  :  pp.  165-168,7]. 
Gray  and  Markel  developed  an  algorithm  which  produces  lattice  realizations  of  filters 
from  the  direct  form  [Ref.  8]. 

The  second  objective  of  the  thesis  is  to  make  the  proposed  ARMA  digital  lattice 
filter  of  [Ref.  4:  p.  662]  adaptive.  An  adaptive  lattice  filter  is  one  in  which  the  lattice 
coefficients  are  automatically  adjusted  by  an  adaptive  algorithm  to  yield  the  optimum 
filter  design.  The  adaptive  lattice  algorithm  derived  in  this  thesis  is  based  on  the  widely 
used  least  mean  square  (LMS)  algorithm.  Adaptive  filters  have  many  applications  [Ref. 
9:  pp.  7-31]  including. 

1.  System  identification. 

2.  Digital  representation  of  speech. 

3.  Adaptive  auotoregressive  spectrum  analysis. 

4.  Adaptive  detection  of  a  signal  in  noise  of  unknown  statistics. 

5.  Echo  cancellation. 

6.  Adaptive  line  cnhancment. 

7.  Adaptive  beamforming. 

The  need  for  an  adaptive  filter  is  made  apparent  by  considering  a  filter  of  fixed  design 
which  is  optimized  for  given  input  conditions.  In  practice,  the  complete  range  of  input 
conditions  may  not  be  known  or  could  change  from  time  to  time.  A  filter  of  fixed  design 
would  not  produce  optimum  results  under  these  conditions.  An  adaptive  filter,  which 
yields  optimum  results  given  changing  input  conditions,  will  give  superior  performance 
to  one  of  fixed  design. 

The  last  objective  is  to  analyze  convergence  properties  of  the  derived  adaptive  lattice 
algorithm.  This  is  accomplished  by  computer  implementation  of  the  adaptive  algorithm. 
The  output  of  a  known  transfer  function  is  compared  to  the  output  of  the  adaptive 
lattice  filter  given  a  common  input.  Plots  of  the  error  between  the  two  outputs  and 
lattice  coefficient  convergence  are  obtained. 

B.  ORGANIZATION  OF  THE  THESIS 

Chapter  II  is  designed  to  present  the  ARMA  parameter  estimation  algorithm  and 
ARMA  digital  lattice  filter  proposed  in  [Ref.  4:  pp.  617-62S].  Computer  simulation  of  the 


algorithm  was  performed  and  results  are  shown.  A  brief  review  of  the  Mullis-Robcrts 
criterion  is  provided  to  establish  a  reference  for  expanding  this  criterion  to  the  proposed 
ARMA  parameter  estimation  algorithm.  An  adaptive  lattice  algorithm  is  derived  in 
Chapter  111  which  makes  the  proposed  ARMA  digital  lattice  filter  adaptive.  The  adap¬ 
tive  lattice  algorithm  is  efficient  and  accurate.  Chapter  IV  contains  experimental  results 
which  show  convergence  aspects  of  the  adaptive  lattice  algorithm  when  applied  to 
ARMA  lattice  filters.  Conclusions  about  the  proposed  ARMA  parameter  estimation 
algorithm  as  well  as  the  derived  adaptive  algorithm  are  discussed. 
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II.  ARMA  DIGITAL  LATTICE  FILTER 

In  this  chapter,  we  will  review  the  Mullis- Roberts  criterion  for  solving  linear  ap¬ 
proximation  problems  and  introduce  analysis  equations  of  the  ARMA  digital  lattice  fil¬ 
ter.  The  criterion  used  in  the  formulation  of  the  ARMA  digital  lattice  filter  is  a 
generalized  form  of  the  Mullis-Roberts  criterion  [Ref.  5:  pp.  227-228J,  which  has  been 
given  as  a  modified  least  mean  square  problem  for  ARMA  parameter  estimation. 

A.  MULLIS-ROBERTS  CRITERION 

The  Mullis-Roberts  criterion  evolved  from  considering  second  order  statistics  in 
conjuction  with  first  order  information  about  a  process  to  obtain  filter  approximations. 
Consider  the  bounded  impulse  response  sequence  h  =  {/?0,  /i,. ...  j  containing  first  order 
information  about  the  filter  h  having  a  frequency  response  function, 


=  y^hke'jojk  (2.1) 

Jfe-0 

Let  {  u,  }  be  a  zero-mean,  unit-variance,  white-noise  sequence  and  { y,  }  be  the  output 
process  corresponding  to  the  input  {  u,  }  ,  then  we  have  the  following  convolutional  re¬ 
lationship  given  by, 


k= 0 


(2.2) 


Second  order  information  about  the  filter  h  is  obtained  from  the  autocorrelation  se¬ 
quence  {  rk  }  given  as 


1 k  Elh}’l+k)  (2-3) 

1=0 

From  equations  (2.1)  and  (2.2),  the  second  order  interpolation  problem  is  to  find  a 
lowest  order  recursive  filter  which  matches  the  data  {  /i„, ... ,  h„.  r0, ... ,  rm  }  ,  where  h,  re¬ 
presents  the  first  order  information  and  r,  the  second  order  statistics. 
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Let  us  now  consider  the  case  where  only  first  order  information  about  a  process  is 
known.  That  is,  given  an  impulse  response  sequence  {  h^,  h„ ...  }  ,  we  want  to  find  a  re¬ 
cursive  filter  of  the  form, 


oc 


Lu 

k= 0 


qp  +  q^  '  +  -  +  qmz 
1  +  axz~x  +  +  anz' 


(2.4) 


which  approximates  H(z)  and  therefore  the  impulse  response  sequence  { /%,  h{, ...  )  .  We 
also  desire  the  frequency  response  to  approximate  the  desired  response  //(e"). 

Suppose  that  //(e")  is  chosen  such  that  it  minimizes  the  integral  squared  error, 


j  m^)  -  =  ||/t  —  A||2  (2.5) 

Using  the  Parseval  relation,  we  can  obtain  an  alternative  definition  of  the  approximation 
error. 


^rJ"|//(O-/VC0)lVar  =  ||/i-/i||2  =S^\hk-hk)2  (2.6) 

k= 0 

If  the  filters  //(;)  and  II(z)  are  driven  by  the  same  white  noise  source,  equation  (2.6) 
describes  the  output  error  between  the  filters  which  we  write  as, 

=  £(y,  -  y,)2  =  £(c,)2  (2.7) 

where  y,  and  y,  are  the  outputs  of  the  respective  filters  when  driven  by  the  same  white 
noise  source  as  in  (2.2).  Minimizing  (2.7)  is  a  nonlinear  programming  problem  requiring 
the  entire  impulse  response  sequence.  As  a  result,  computational  efforts  for  obtaining  a 
solution  arc  inefficient.  A  modification  to  the  problem  was  introduced  [Ref.  5:  pp. 
227-22S]  which  considered  a  cost  function  that  is  quadratic  in  the  coefficients  of  the  re¬ 
cursive  filter  given  by  (2.4).  The  modification  is  described  as  follows.  Let 

m  =  +  •  +  <1rrZ~m)  (2.8) 
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and 


A{z)  =  r  '( 1  +  a,.--1  +  -  +  anz~n)  (2.9) 

be  the  numerator  and  denominator  polynomials,  respectively,  in  (2.4).  where  S—  max 
(nun).  The  task  now  is  to  find  coefficients  which  minimize  the  quadratic  form. 

E{C‘)2  =  17  J  (2.10) 

This  is  a  standard  quadratic  minimization  problem  whose  integral  can  be  expressed  in 
terms  of  the  coefficients  of  polynomials  A(z)  and  Q(z)  and  the  data 
{/io, ... ,  hm.  r0. ... ,  rm},  relating  to  the  filter  //(x).  This  problem  is  shown  in  Figure  1  and 
equation  (2.10)  is  known  as  the  Mullis-Robcrts  criterion. 


Figure  1.  Modified  least  squares  problem 


B.  GENERALIZED  MULLIS-ROBERTS  CRITERION 

In  order  to  define  the  new  criterion  used  for  ARMA  parameter  estimation  we  con¬ 
sider  the  following  transfer  function  with  input  sequence  {  x(k),  k  =  1,2, ...  }  and  output 
sequence  [y(k),  k  =  1,2, ...  }  written  as, 

(2.11) 


H(z)  = 


NyU) 
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where  H}(:)  and  //,(-)  are  reference  polynomials  which  we  desire  to  model.  An  equivalent 
model  is  shown  in  Figure  2  where  u(k)  is  an  intermediate  signal,  and  the  realization  is 
similar  to  that  of  the  direct  form  realization  II  [Ref.  10:  p.  151].  Let  the  intermediate 
sequence  {  i<{k),  A  =  1.2. ...  }  be  a  zero-mean  white  gaussian  process.  The  model  of  Fig¬ 
ure  2  can  then  be  transformed  into  the  model  of  Figure  3  with  u(k)  as  the  common  input 


Figure  3.  Transformed  model 


Note  that  we  earlier  defined  transfer  functions  H,{z)  and  HM{z)  as  polynomials  of  the 
reference  model.  For  each  transfer  function  an  estimation  polynomial  is  defined  such 
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that 


/1(c)  =  1  +  UjZ  '  + 

,  —5 

‘  +  GfZ 

estimates  JIx(z) 

(2.12) 

7?(r)  =  +  b\Z  1  +  ■ 

•+V-" 

estimates  //v(z) 

(2. 13) 

where  a,  and  b„  (i  =  1,2 . s)  and  (J  =  1,2. ... ,  r)  ,  are  the  AR  and  MA  parameters,  re¬ 

spectively.  of  the  combined  ARMA  model  formed  by  A(z)  and  B(,z).  This  refined  least 
squares  problem  is  shown  in  Figure  4  and  is  a  generalized  form  of  the  Mullis- Roberts 
criterion.  If  the  reference  model  polynomial  !!,(:)  in  Figure  4  is  equal  to  unity,  we  obtain 
the  Mullis-Roberts  criterion  shown  in  Figure  1.  Therefore,  by  including  reference 
polynomial  7/,(r),  the  new  criterion  for  ARMA  parameter  estimation  becomes, 


Es,  =  -^-i\Hv(z)A(z)- Hx(z)D(:)\2  4-  (2.14) 

+.7ZJ  J 


u(k) 

y(k) 

Hy(Z) 

A(z) 

+i  ®(k> 
©— * 

Hx<2> 

x(k) 

B(z) 

•f 

Figure  4.  Refined  least  squares  problem 


Minimizing  EtJ  is  accomplished  by  calculating  the  coefficients  of  A(z)  and  B(z)  which 
minimize  e(ky  of  Figure  4.  Another  form  of  eq  (2.14)  is  obtained  by  applying  Parseval's 
theorem  and  is  expressed  as, 

=  El(A(zlv(k)  -  Z?(r).v(*))2]  (2.15) 

which  is  obvious  from  Figure  4.  The  coefficients  a„  (i  =  1 . s),  and  br  {j  -  1, ...  ,  t), 

which  minimize  can  then  be  calculated  using  the  normal  equations  for  ARMA 


8 


parameter  estimation.  In  order  to  obtain  the  normal  equations  for  the  problem  in  (2.15), 
let  us  define  the  following: 


av  =  O)  ...  aj  and  bv=  [ft,  ...  b,}  (2.16) 

are  the  vectors  of  AR  and  MA  parameters,  respectively, 

M*)  =  l  ><*)  -y(k  -  S)  —x(k) ...  -x(k  -  i)  ]  (2.17) 

is  the  data  vector  consisting  of  both  input  and  output  data  elements  and 

RSit=  E[hs/kfhjk)]  (2.18) 

is  the  data  autocorrelation  matrix.  The  criterion  is  to  minimize  the  mean  squared  error 

£[V(A')  ]-£[[[  1  K,]!] 

L  (■>  i9) 

=  £[[><*)+  C  aj,f  *o  K'  ]  ]2] 

Now  following  the  standard  calculus  of  variation  optimization  procedure  for  minimizing 
£v,  (Ref.  11:  pp.  100-110],  yields  the  normal  equations  in  the  matrix  form, 

[  1  aw  b0  b,f  ]  Rjr  =  [  min  Est  0  0  0  ]  (2.20) 

It  is  interesting  to  note  that  if  £„  in  equation  (2.14)  is  zero,  then  we  have  the  following 
equality, 


H{z)  = 


//v(-) 


Biz) 

A(z) 


(2.21) 


so  the  estimate  for  the  total  reference  model  H(z)  is  the  ratio  of  B(z),  the  MA  part  and 
A(z),  the  AR  part  of  an  estimated  ARMA  model. 

C.  ARMA  PARAMETER  ESTIMATION 

Now  that  the  criterion  for  ARMA  parameter  estimation  has  been  established,  the 
solution  method  to  estimate  the  ARMA  model  parameters  to  minimize  equation  (2.14) 
is  considered.  Let  x(k)  and  y(k)  be  the  input  and  output  signals,  respectively,  of  the  es¬ 
timated  ARMA  model.  Using  a  difference  equation  representation,  this  process  is  de¬ 
scribed  by 


(2.22) 


S  I 

y(k)  =  -  y  a7-  y\k  -j)  +  y  6,  x(k  -  i ) 

7=1  (=0 

For  these  input  and  output  signals  we  define  four  estimation,  or  prediction,  models  as 
follows.  The  forward  estimation  signal  for  x{k)  is  defined  as, 

t  s 

x/k)  =  ~YJbi  *(k  ~  0  +  X  aJ  ><*  (2'23) 

/=i  j=\ 

where  b;  and  af  are  the  corresponding  estimation  parameters.  The  forward  esitmation 
signal  fory(A)  is  similarly  defined  as, 

S  t 

y/.k)  y{k  -j)  +  Yj  b‘  x(k  ~  '■)  (2-24) 

7=1  i=l 

The  backward  estimation  errors  for  x(k  —  /)  and  y(k  —  s)  are  then  given  by, 


/-I  5-1 

.ib(k -t)  =  ~Yb?  x(k -i)  +  Y af  -J) 

(=0  7=0 


(2.25) 


5-1  [-1 

Sb(.k  -  s)  =  -  y  af  y{k  -j)  +  y  bf  x(k  -  /) 


7=0 


i=o 


(2-26) 


where  the  superscripts  g  and  d  indicate  the  backward  estimation  parameters  for  .r  and 
y.  respectively. 

From  these  estimation  models,  we  can  now  obtain  the  four  prediction  errors  at  any 
given  time  k.  These  errors  are  expressed  as  differences  between  the  predicted  value  and 
actual  value  of  an  input  or  output  signal,  namely, 


vUk)  =  -x{k)  +  Xj{k) 

(2-27) 

<,{k)=y{k)~yM) 

(2.28) 

gs,t(k)  =  x(k  xb{k  -  t ) 

(2.29) 

1 


I 
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ds/k)  =yik  -  s)  -}'b(k  -  s) 


(2.30) 


We  now  use  the  vector  notation  to  simplify  the  expressions  for  prediction  errors.  In  the 
following,  the  forward  error  elements  corresponding  to  both  jr  and  y  form  a  vector 
vJA),  given  by 

VW  =  [  -v*,(k)  =  <r  (2.31) 

and  the  backward  error  vectors  are  given  by 

gSJ(k)  = -hs/k)Gl  (2.32) 

djk)  =  hJA)  d£  (2.33) 

where  Cti,  and  G,,r  and  D  ,  arc  the  forward  estimation  parameter  matrix  and  backward 
estimation  parameter  vectors,  respectively,  defined  as 


'0 

X  * 

a>  ...  a , 

1  b f 

... 

Cw- 

1 

1  5 

a]  ...  a] 

0  b\ 

...  b>_ 

(2-34) 

GJ./  =  [4 

...  flf-i  0 

bi  - 

bU 

1] 

(2.35) 

c 

1 — 1 

II 

O 

i 

)  • 

d  , 

-  1 

bo  - 

bU 

0] 

(2.36) 

It  can  be  shown  that  the  prediction  errors  satisfy  some  orthogonal  conditions.  These 
conditions  are  similar  to  those  found  in  AR  modeling  problems  [Ref.  6  :  pp.  116-121]. 
We  now  list  the  orthogonality  conditions  for  the  ARMA  formulation  in  discussion 
without  proof  as  the  following: 

E[  <,(*)  y(k  ~j)  ]  =  0  £[  ty»  y(k -j)  ]=0 
£[  <f(£)  x(k  -  0  ]-0  £[  ii,(k)  x(k-i)  ]  =  0 

Eigs/k-U  y(k-jll  =  0  (3.37) 

Etgjk  -  1)  x(k  —  /)3  =  0 
Eldjk-l)  y(k  -7)]  =  0 
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£(X.,(*-1)  x{k  —  /)]  =  0 


where  i  =  1,2, ... ,  /  and  j  =  1,2, ... ,  s. 

In  Section  B  we  have  obtained  a  set  of  normal  equations  in  terms  of  R, ,  and  the 
ARMA  estimation  parameters.  In  this  section,  we  have  defined  four  sets  of  forward  and 
backward  estimation  parameters  and  established  some  orthogonal  relationships.  In  what 
follows,  we  derive  a  set  of  equations  which  relates  the  coefficients  of  the  estimated 
AR.MA  model  with  those  of  the  forward  estimation  parameter  matrix  C,,.  Consider  the 
expected  value  of  the  forward  prediction  error,  vI  ((A)  and  the  data  h,,,(/c).  Since  the  pre¬ 
diction  error  is  orthogonal  to  all  past  samples  of  data  y(k  —  j),  x(k  —  i)  but  not  to  y(k) 
or  x(k)  as  listed  in  (2.37),  the  result  is  a  matrix  which  is  defined  as 

£[v(*)  M*)]=  .  :  n  .  n  (2.38) 

J  L  0  Cj  0 

where 

*  1  =  ~E  [  <,(*)>(*)  ]  =  -E  [  i&k)  rl(k)  ]  =  Eff  (2.39) 

is  the  crosscorrclation  between  the  forward  prediction  errors  of  x  and  y  at  lag  zero. 

c2  =  £  [  v*,(k)  x(k)  ]  =  £  [  (r*,(A))2  ]  =  (2-40) 

is  the  forward  prediction  error  power  of  x{k) 

=  E[^,(k)y(k)  ]  =  £  [  «,(*))2  ]  -  £j,(  (2.41) 

is  the  forward  prediction  error  power  of  j(A)  and 

-  -£  [  X(k)  >  -£  [  tijik)  \,{k)  ]  =  EfJ  (2.42) 

is  the  crosscorrelation  between  the  forward  prediction  errors  of>'  and  x  at  lag  zero.  In 
another  interpretation,  the  left  hand  side  of  equation  (2.38)  can  be  written  as, 

£  [  xSJ{k)T \,{k)  ]  =  £  [  C„  hjkfhjk)  ]  =  CIf  Rsl  (2.43) 

and  we  have 
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Cj, i  Rj,r  ~ 


Exy 

‘-s.l 

Ey 

s.t 


Ex 

us.t 


-i,r 


(2.44; 


from  equation  (2.38). 

A  similar  approach  for  both  backward  prediction  errors  yields  the  following 


0 

Egd 

~ S,l 

0 

E§ 

0 

Z’SJ 

0 

Egd 

(2.45) 


In  order  to  express  the  coefficients  of  the  ARMA  estimation  mode)  in  terms  of  the  co¬ 
efficients  of  the  parameter  estimation  matrix  Cs!  and  parameter  vectors  G, ,  and  D!-(  , 
consider  the  combination  of  the  normal  equation  (2.20)  and  the  parameter  estimation 
matrix  and  vectors.  From  equation  (2.44)  the  normal  equation  for  ARMA  parameter 
estimation  may  be  written  as 


'o 

as.t 

1 

*£/ 

R, ,  = 

* 

0 

Ex 

s.t 

0  ’ 

1 

0 

blr 

* 5,1 

0 

Exy 

u s.t 

0 

(2.46) 


Combining  equation  (2.46)  with  (2.20)  we  obtain 


’  1 

a,./ 

*0 

Z^min 

0 

0 

o' 

0 

a;. 

1 

K, 

R,  = 

Ex> 

0 

Ex 

‘'S.t 

0 

1 

aL 

0 

Ey 

s.t 

0 

Exy 

C'S.l 

0 

In  equation  (2.47)  multiply  row  2  by  then  subtract  row  2  from  row  3  and  equate 
the  resuit  to  row  1.  We  then  obtain  the  following  relations 


Jj,/ —  a5,t  &S.t 


Ex 


b0  - 


Exy 

Fx 

‘-'S.t 


xy 


E  J 

£s, 


(2.4S) 


(2.49) 


(2.50) 


and, 
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Est  min  =  , 


■xy.2 


E 

*-*5.1 


(2.51) 


Calculation  of  the  ARMA  estimation  model  coefficients  requires  knowledge  of  the  esti¬ 
mation  parameters  a*„  a>,„  b;„  tyr,  b0  and  the  values  £*,,  E>„  E'j  .  Recursive  formulas 
calculate  the  estimation  parameters  as  the  AR  or  MA  order  of  the  estimation  model  in¬ 
creases  by  one.  Given  the  parameters  of  the  ARMA  estimation  model  of  order  (s ,  /) 
these  formulas  calculate  the  (s+  1  ,  r)  or  (s,  t+  1)  ARMA  parameters.  For  a  compre¬ 
hensive  derivation  of  these  recursive  formulae  see  [Ref.  4:  pp.  6 1 9-62 1 J.  The  AR-type 
recursive  formula  for  the  forward  parameter  estimation  matrix  and  backward  estimation 
parameter  vectors  is 

Gi+l,r  ~  Cj  j  Ij  +  U]  Dj  ,  12 

G*+l,r  -  [  GJ./  +  u2  Dj,,  ]  *1  (2-52) 

=  ^2  "h  d  U3  C 5  l  GJ  (  4-  U5  Di  f  ]  I  j 


where 


».--(4r'u  o 

_  p?  4 


u2=- 


EJ 

U3  =  -[r,  r j]£T; 

(Elf  r4  -  £?,  t3) 


= 


u5  =  u4  «2 


and  the  (s+  l,t)  prediction  error  powers  are  recursively  calulated  as  follows 

Es+\,t=  Esj  +  UI C  t,  r2  ] 

El+\.,  “  +  u2  *4 

£?+],,  -  Ef,t  +  C  T,  x2  Hu[  +  u4r3  +  u5 r4 

£?+u  =  Eh  +  m2 


(2.53) 


(2.54) 


The  matrices  I,  and  I2  are  of  dimension  (5  4-  t  +  2  x  j  4-  t  4-  3).  They  are  introduced  to 
provide  symmetry  to  the  matrix  algebra  and  preserve  initial  condition  calculations,  We 
design  the  matrices  to  perform  the  following  operations 


14 


C  a>! 

...  a>J+1  o)s+2 ...  tuI+f+2  ]  I; 

II 

1 — 1 

S 

•  w*+l  0  “h+2 

■  coj+i+2  3 

(2.55) 

1  u\  • 

■  ^h+l  W4+2  •••  wi+/+2  3 

II 

1 — 1 

0 

e 

•“Vt-l  0(0 s+ 2- 

•  Wj+H-1  ] 

(2.56) 

The  values  t,  through  t4  can  be  obtained  through  the  correlation  data  and  forward  and 
backward  estimation  parameters.  We  express  them  mathematically  as 

L  t,  r 2  D  =  £  [  ds  t{k  -  1)  \sl(k)  ] 

r3  =  ~E  C  y(k  -  s  -  1)  gsl(k)  ]  (2.57) 

t4  =  EZyik-s-  \)dsJ,k)  ] 

The  MA-tvpe  recursive  formula  for  the  forward  estimation  parameter  matrix  and  back¬ 
ward  estimation  parameter  vectors  is  obtained  in  a  similar  manner.  The  recursive  for¬ 
mula  is  given  by 

Ci,r+1  =  c*.r  I3  +  nT Gs,t  I4 

Dv+1  =  [D^  +  «2Gw]l3  (2.58) 

Gj,r+1  =  Gs,/  *4  +  [  n3  Cs,t  +  rt4  D1>(  +  «5  Gs  t  ]  I3 


where 


“1  =  ~  ’  C  t'2  ] 

«2  = 


r"d 

L'*,t 


Eh 

n3  =  ~  [  t'2  ]  E -j 

(Ef/T'4-E*fT'3) 

*  [  {Efff  -  El,  ] 

«5  =  «4  n2 


(2.59) 


and  the  (s,t+  I)  prediction  error  powers  are  calculated  using  the  following  recursive  for¬ 
mulas 


—  ESt,  +  flj  E  T  1  T  2  J 
Ef,t+\  ~  r>3  +  ni  T'a 

El,+i  =  El,  +  [  t',  t'2  ]  n3r  +  «4  t'3  -l-  «5  r'd 


ud 


(2.60) 


where  I3  and  I4  are  (5  +  /  +  2  x  s  +  1  +  3)  dimensional  matrices  which  we  have  designed 
to  perform  the  following  operations 


C  Wj  .. 

'  1  Us+2  ...  Ws+(+2  H  I3  —  C  w!  •••  ^J+l  ••• 

<»s+t+2  0  J 

(2.61) 

C  O)!  . 

-  <uJ+,  cos+2  ...  cos+t+2  J  I4  =  [  0  O),  ...  Wj  0  0)^2 

...  J 

(2.62) 

The  values  of  t',  through  t'4  are  calculated  using  correlation  data  in  conjunction  with 
current  forward  and  backward  parameter  estimation  values.  We  express  these  quantities 
mathematically  as 

[t',  t'2  J  =  —E  [  gs  t{k  —  1)  vx  f(A)  ] 

r'3  =  -£[*(*-/-  1K./A)]  (2.63) 

T\  =  El.x(k-t-  1  )£,.,(*)] 

It  is  interesting  to  note  that  the  .MA-type  recursive  formula  is  the  complimentary  form 
of  the  AR-type  formula  and  that  the  two  are  identical  if  the  variables  associated  with  the 
input  signal  .v(/;)  and  the  variables  associated  with  the  ouput  signal y(k)  are  interchanged. 
That  is,  we  replace  vv(»c),  Gw  and  g,-r(A)  with  -x(k).  D;>,  and  ds,(k)  and  vice  versa. 

1.  Experimental  Results 

The  ARM  A  parameter  estimation  algorithm  of  [Ref.  4:  pp.  6 1 9-62 1]  based  on 
the  recursive  formulas  of  equations  (2.52)  and  (2.58)  was  implemented  using  the  Fortran 
program  found  in  Appendix  A.  This  program  calls  subroutines  which  compute  the 
ARMA  model  parameters  as  the  AR  order  is  increased  by  one  and  as  the  MA  order  is 
increased  by  one.  These  subroutines  are  shown  in  Appendix  B  and  Appendix  C.  respec¬ 
tively.  In  the  main  program,  an  input  data  sequence  of  white  Gaussian  noise  is  passed 
through  a  known  reference  model  producing  an  ouput  data  sequence.  We  obtain 
autocorrelation  and  crosscorrelation  data  from  these  input  and  output  sequences.  The 
correlation  data  is  used  to  calculate  initial  values  of  the  error  powers  for  x  and  y  as  well 
as  t,  through  t4  and  t',  through  r'4 .  Next  we  obtain  estimates  of  the  reference  model  by 
employing  the  recursive  formulas  (2.48)  through  (2.50),  (2.52)  and  (2.58).  Several  refer¬ 
ence  models  were  estimated  beginning  with  a  strictly  AR  process  of  order  s  =  4  having 
as  its  transfer  function 


m 


_ 1 _ 

1  —  0.2  z_l  +  0.62  z~2  -0.152  r-3  +  0.3016  z-4 


(2.64) 
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j 


The  actual  values  of  the  reference  model  parameters  and  the  ARMA  model  parameters 
which  estimate  this  reference  model  are  listed  below 


ACTUAL 

ESTIMATED 

0.2000 

0.2005831 

-0.6200 

-0.6207655 

0.1520 

0.1527565 

-0.3016 

-0.3020376 

1.0000 

1.0001506 

We  next  consider  a  second  reference  model  with  MA  order  /  =  2  and  AR  order  s  =  3 
havina  transfer  function 


//U-)  = 


0.5  —  0.40  c-1  +  0.S9  r~2 
1  -  0.20  2~x  -  0.25  r-2  +  0.05 


(2.65) 


The  true  reference  model  parameters  and  ARMA  model  parameter  estimates  are  shown 
to  be 


AR: 


MA: 


ACTUAL 

ESTIMATED 

0.2000 

0.1993060 

0.2500 

0.2496567 

-0.0500 

-0.0491961 

0.5000 

0.5002602 

-0.4000 

-0.3997071 

0.8900 

0.SS94749 

A  third  example  with  MA  order  t  =  2  and  AR  order  5  =  4  having  transfer  function 


H(z) 


_ 1  4-  0,2  r~'  —  0.99  r~2 _ 

I  —  0.2  r-i  +  0.62  z~2  —  0.152  z~3  +  0.3016  r_J 


(2.66) 
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was  considered  for  which  we  obtained  the  following  actual  and  estimated  reference 
model  parameter  values 


ACTUAL 

ESTIMATED 

AR:  0.2000 

0.2011805 

-0.6200 

-0.6223S03 

0.1520 

0.1534197 

-0.3016 

-0.3036S23 

MA:  1.0000 

0.999763S 

0.2000 

0.1998342 

-0.9900 

-0.9SS6S52 

We  consider  as  a  final  example  the  reference  model 
i  =  3.  respectively,  with  specific  transfer  function 

//(*)  = 

0.5  —  0.95  z-1  +  1.33  2 

1  +  1.69  z"1  -0.962 

The  actual  and  estimated  ARM  A  parameters  are 

ACTUAL 

ESTIMATED 

AR:  -1.6900 

-1.6981325 

0.9620 

0.969099S 

-0.2000 

-0.2018440 

MA:  0.5000 

0.4995653 

-0.9500 

-0.9553509 

1.3300 

1.3346767 

-0.9790 

-0.9864898 

of  AR  order  and  MA  order  5=3  and 


~2  -0.979  z~3 

z_:  +  0.2  z~3 


(2.67) 


The  above  examples  demonstrate  the  validity  of  the  parameter  estimation  algo¬ 
rithm  of  (Ref.  4:  pp.  619-621 J.  Many  reference  models  were  estimated  using  this  algo¬ 
rithm.  including  pure  MA  processes,  for  which  accurate  estimates  were  obtained. 
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D.  LATTICE  STRUCTURE 

In  section  C  \vc  developed  expressions  for  the  forward  and  backward  prediction  er¬ 
rors.  namely,  those  of  equations  (2.31),  (2.32)  and  (2.33).  From  these  prediction  error 
equations  we  can  design  elementary  AR,  MA  and  AR.V1A  lattice  structures  or  sections. 
Each  elementary  section  satisfies  the  orthogonal  conditions  as  listed  in  equation  (2.37). 
From  the  prediction  error  recursive  formulae,  equations  (2.31).  (2.32)  and  (2.33),  we 
construct  the  AR-type  elementary  lattice  section  as  follows.  Consider  the  following  data 
set  of  order  (s  +  1 ,/)  consisting  of  input  and  output  data  elements, 

W*)  If  =  CM*)  ...y(k  -  s )  -x(k) ...  -x(k  -t+  1)  -x(k  ~  ')  ]  (2.68) 

W*)  if  =  CM*  -  1)  ...y(k  -s-  1)  -x(k  -  1) ...  -x(k  -  t)  0  ]  (2.69) 

where  If  and  If  are  the  transposes  of  the  matrices  I,  and  I2  defined  in  equations  (2.55) 
and  (2.56).  We  obtain  a  recursive  relationship  between  the  forward  prediction  errors 
v(A)  of  order  (s-r  I.r)  and  order  ( s,t )  by  substituting  equations  (2.52).  (2. 68)  and  (2.69) 
in  equation  (2.31 )  such  that 

'*+ 1./^)  =  I*i+i. f(^)  Cj+i.i 

=  h„,  ,«)[irC,r,  +  I;D,r,u,  ]  (2.70) 

=  V, ,,{k)  ~  4  ,(A-  -  l)u, 

The  backward  prediction  error  recursions  are  obtained  in  a  similar  manner  and  the 
AR-type  error  recursions  are 

wM)  =  -  u\  dsJk  ~  i) 

v j+ 1 M •'  =  ls.r  0 k )  +  t{\  dsAk  ~  I)  ( ?  j 

g:+<Jk)  =  gSJ(k)  -  u-,  dS  I(k) 

4+  !,:(*)  =  dsAk  “  I)  +  [  ui  ul  ]  ys/kf-  u-i  Ss+\M) 

where  u,  =  C  u\  n]  ]  and  u3  =  C  u)  f3-  The  AR-type  elementary  lattice  inverse  section 
based  on  these  error  recursions  is  shown  in  Figure  5. 


! 
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Figure  5.  AR-type  elementary  lattice  inverse  section 

We  design  the  MA-typc  elementary  lattice  inverse  section  in  a  similar  manner  using  the 
following  representation  of  the  data  set  of  order  (s.i  +  1) 

h,  ,+  l  (A)  lj=[  y(k)  y(k  -  1) ...  y(k  -  s)  ~x(k) ...  ~x(k  -t- f  1)  -x(k  -  t)  ] 
h„+)  (A)  lj=  [  y(k  -  1)  ...y{k  -  s)  0  -x(k  -  1) ...  -,r(A  -  /  -  1)  ] 

and  by  substituting  equations  (2.5S)  and  (2.12)  into  equation  (2.31).  we  obtain  the  for¬ 
ward  prediction  error  recursion  as  the  MA  order  is  increased  by  one,  namely, 

'*/+  ,(A)  =  v*,(A)  +  gs/k  -  1)  7 , 

=  lir  W  -  wl  SsJA  ~ 

The  backward  prediction  error  relationships  are  obtained  in  a  similar  manner  and  are 
given  by 

^s,t+  i(A)  =  ds,{k)  /?2  gj./A) 

StJ*  i(A)  =  Ss/k  -  1)  -  [  n*  Wj  ]  \s/k)T-  «4  dsl+](k)  *  ‘ 

The  MA-type  elementary  lattice  inverse  section  based  on  these  error  recursions  is  shown 
in  Figure  6. 
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We  now  construct  the  ARM  A  elementary  lattice  section  from  the  AR  and  MA 
prediction  error  recursions.  Assuming  that  the  prediction  errors  are  known  for  a  given 
model  order  (s.t).  the  (s-r  l.t)  prediction  errors  can  be  calculated.  These  prediction  errors 
of  order  (s  +  l.t)  are  then  updated  as  the  MA  order  increases  by  one  resulting  in  a  pre¬ 
diction  error  of  order  (s-r  l.t  +  1  i.  We  consider  the  forward  prediction  error  for  x  as  the 
AR  order  is  increased  by  one.  specifically 

1) 

Now.  ,  (A)  becomes  the  current  value  of  the  forward  prediction  error  for  .v  and  when 
we  calculate  the  (s  +  l.t  +  1 )  forward  prediction  error  we  have  from  equation  (2.73) 

h+u+i  (A)  -  &i  j  (*)  +  ”i  Ss-ru  (A  ~  1)  (2"6) 

Equation  (2.76)  can  be  expressed  in  terms  of  the  (s.t)  forward  prediction  errors  of  x  by 
making  appropriate  substitutions  for  v;_,,  (A)  and  g„ U(A  —  1).  That  is,  we  substitute 
v;.Ul  (A)  and  g,_u(A  -  1)  of  equation  (2.71)  in  equation  (2.76)  to  obtain  the  (s+  l.t+  1) 
forward  prediction  errors  for  .r,  namely. 


W,  (k)  =  vs,i  (A)  -  uUu  (A  -  0  +  « i  [  (A  -  1)  -  u2  d5J  (A  -  1)  ]  (2.77) 


Grouping  the  terms  we  obtain 


'•i+i,,+i  (A)  =  vf,  ( A )  -  (uf  +  «f  u2)  ds  ,  (A  —  1)  +  rtf  gs  t  (A  -  1) 


(2-78) 


The  forward  prediction  error  recursion  for  y  of  order  (s  +  l,t+  1)  is  obtained  in  a  similar 
manner.  We  begin  with  the  (s+  l,t)  order  update  of  the  prediction  error  and  after  it  is 
computed,  update  the  MA  order.  Specifically,  we  have 

=  +  ^  <,(*-!)  (2-79) 


and  from  equation  (2.73) 

1  (*)  =  tf+u  (k)  -  «f  &+lfl(A  -  1)  (2.S0) 

Substituting  (2.71)  for  v;.,  t  (k)  and  g,.h,(k—  1)  in  equation  (2.80)  then  grouping  terms 
we  obtain  the  ($  +  l,t  +  1)  forward  prediction  error  recursion  for y, 

vf+i,/+i  (A)  =  4,i  (A)  +  (4  +  “ 2 )  dsj  (A  ~  1)  ~  «?  Ss,,  (A  ~  1)  (2-81) 

The  (s+  l,t+  1)  backward  prediction  errors  for  x  and  y  are  derived  in  a  similar  manner 
and  are  given  by, 

&+U+1  (A)  =  Ssj  (A  ~  1)  +  («3  +  ui)  vsj  (*)  ~  («3  +  «4  *4)  (k) 

ds+u+ ,  (k)  =  <,  (A  -  1)  -  uf  vf,  (A)  +  uf  4,  (A)  (~  ,fc) 

The  ARMA  elementary  lattice  inverse  section  is  shown  in  Figure  7  where  the  coefficients 
are  related  to  the  prediction  error  recursions  by  the  following 

wj  =  (uf  +  rtf  u:),  wj  =  rtf,  vvj  =  (tq  +  rtf  u2),  ivf  =  rtf  (2.83) 

*5  =  (rtf  +  «4  uf),  wl  =  (rtf  +  «4  uf),  uf  =  uf,  Wg=uf  (2.84) 


X  X 


Figure  7.  ARM  A  elementary  lattice  section 

Wc  see  from  Figure  7  that  each  elementary  ARMA  lattice  inverse  section  contains  eight 
coefficients. 

From  the  AR,  MA  and  ARMA  elementary  lattice  inverse  sections,  we  can  obtain 
synthesis  lattice  structures.  These  structures  provide  a  means  of  working  with  lattice  re¬ 
alizations  as  linear  filters.  The  resulting  AR,  MA  and  ARMA  elementary  synthesis 
lattice  filters  are  shown  in  Figures  S  and  9  respectively. 
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Figure  8..  Top 


t 


> 


Figure  9.  ARMA  elementary  lattice  filter 


Summarizing,  in  this  chapter  we  have  reviewed  the  Mullis- Roberts  criterion,  intro¬ 
duced  the  ARMA  parameter  estimation  as  a  generalized  Mullis-Robcrts  criterion  and 
obtained  analysis  and  synthesis  forms  oflattice  structures.  We  notice  that  each  ARMA 
elementary  lattice  section  consists  of  eight  reflection  parameters  and  the  calculation  of 
these  parameters  requires  the  autocorrelation  and  crosscorrelation  information  as  ob¬ 
tained  from  the  input  output  data  of  the  reference  model.  Also,  we  obtained  a  set  of 
equations  relating  the  final  model  estimation  parameters  and  the  prediction  error  model 
parameters  which  inturn  are  obtained  from  the  eight  lattice  parameters. 
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III.  ADAPTIVE  LATTICE  ALGORITHM 

A.  LEAST  MEAN  SQUARE  ALGORITHM 

The  study  and  design  of  adaptive  filters  is  known  to  be  a  very  important  part  of 
statistical  signal  processing.  Many  adaptive  algorithms  have  been  developed  to  support 
the  application  of  adaptive  filtering  in  communications  and  control  [Ref.  12J.  An  adap¬ 
tive  filter  is  characterized  by  the  ability  of  its  filter  coefficients  to  adjust  (self-optimize) 
automatically  and  yield  an  optimum  filter  design.  Two  processes  occur  within  an  adap¬ 
tive  filter,  namely,  the  adaptation  and  the  filtering  processes.  During  the  filtering  process 
a  desired  signal  is  applied  to  an  adaptive  algorithm  as  a  reference  for  adjusting  the  filter 
coefficients.  Figure  10  shows  a  block  diagram  of  the  adaptive  modeling  process.  Refer¬ 
ring  to  Figure  10,  let  >•(/:)  be  the  output  of  the  filter  at  time  k.  By  comparing  the  output 
with  the  desired  signal  d(k)  ,  an  error  signal  e{k)  is  generated.  The  adaptive  algorithm 
of  the  filter  uses  this  error  signal  to  generate  corrections  which  are  applied  to  the  filter 
coefficients  such  that  an  optimum  solution  is  obtained.  An  optimization  technique  called 
the  method  of  steepest  descent  provides  an  approach  to  solving  this  problem.  The  pro¬ 
cedure  is  as  follows: 

1.  Assign  initial  values  to  all  filter  coefficients. 

2.  Using  these  initial  values,  compute  the  gradient  vector,  whose  individual  elements 
equal  the  first  derivatives  of  the  mean-squared  error  with  respect  to  the  filter  coef¬ 
ficients. 

3.  Compute  values  for  the  filter  coefficients  by  changing  the  initial  values  in  the  di¬ 
rection  opposite  that  of  the  gradient. 

4.  Return  to  step  2  and  repeat  the  procedure. 

There  is,  however,  a  limitation  to  this  procedure.  The  steepest  descent  algorithm  requires 
exact  measurements  of  the  gradient  vector  at  each  iteration  which,  in  practice,  is  not 
possible.  Therefore,  the  gradient  vector  must  be  estimated  and  consequently,  errors  are 
introduced.  An  algorithm  is  required  which  computes  the  gradient  from  the  available 
data.  The  least  mean  square  (LMS)  algorithm,  developed  by  Widrow  and  HofT,  is  widely 
used  and  is  very  convenient  to  implement  in  real  time  hardware  [Ref.  13:  pp.  96-104]. 
Let  y{k)  be  the  output  of  the  filter  and  d(k)  the  desired  signal  at  time  k  as  shown  in 
Figure  10.  We  compute  the  error  by  taking  the  difference  between  these  two  signals, 
namely, 
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Figure  10.  Adaptive  modeling  block  diagram 


c(k)  =  J(k)  -  y(k) 


(3.1) 


The  value  of  the  mean-squared  error  is  the  expected  value  of  the  error  squared, 
e-(k)  ]  and  the  gradient  vector,  V  (k)  ,  is  the  first  derivative  of  the  mean-squared  er¬ 
ror.  The  gradient  vector  is  given  by 


CM(/C) 


e(k) 


(3.2) 


where  w (k)  is  the  time  dependent  filter  coefficient  vector.  The  recursion  for  the  filter 
coefficient  which  changes  the  old  value  in  the  direction  opposite  to  that  of  the  gradient 
is  then  given  by, 

„(*)_„(*_  d+X  „[  -v(t)] 


where  >\(k)  is  the  filter  coefficient  vector  estimate  at  the  k rt  iteration,  v,(k  —  1)  is  the  past 
filter  coefficient  vector  estimate,  n  is  the  convergence  (gain)  constant,  e(k)  is  the  error 
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signal  at  the  k,h  iteration  and 


c\s(k) 


e{k)  is  the  instantaneous  gradient.  The  implemen¬ 


tation  of  this  algorithm  proceeds  as  follows: 

1.  Assign  initial  values  to  the  filter  coefficients. 


2.  Compute  the  value  for  the  error  signal  e(k). 

3.  Calculate  the  updated  estimate  of  the  filter  coefficients  using  the  instantaneous 
gradient. 

4.  Increment  the  time  index  by  one  and  return  to  step  1. 


Convergence  properties  of  the  LMS  algorithm  are  well  documented  within  the  literature. 
The  choice  of  a  gain  constant  /i  is  arbitrary  however,  theoretical  bounds  have  been  de¬ 
rived  for  n  ,  given  by  [Ref.  14:  pp.  101-106], 


0  <n  < 


2  1 
'max  '  TrCR^] 


(3.4) 


where  /m„  is  the  maximum  eigenvalue  of  the  input  autocorrelation  matrix,  R„,  and 
where  Tr  [  R„  ]  is  the  trace  of  the  matrix  R„. 

B.  DERIVATION  OF  THE  ADAPTIVE  LATTICE  ALGORITHM 

The  adaptive  lattice  algorithm  developed  in  this  thesis  uses  concepts  of  the  LMS 
algorithm  discussed  in  section  A  and  applies  them  to  the  ARMA  digital  lattice  filter 
proposed  in  Chapter  II.  Consider  the  ARMA  digital  lattice  filter  of  Figure  11,  which 
consists  of  two  cascaded  elementary  lattice  sections.  The  filter  coefficients  (weights)  are 
defined  such  that  ivf  represents  the  i"'  lattice  coefficient  at  stage  m  of  the  lattice  structure. 
In  this  figure  we  have  a  two  stage  lattice  and  there  are  eight  coefficients  per  elementary 
lattice  section.  The  output,  y(k),  of  the  lattice  filter  can  be  determined  from 


=  e)\(k)  +  ivj  e^(k  -  1)  -  *•’  eljk  -  1) 


Forward  errors  at  a  given  stage  m  of  the  lattice  filter  are  defined  as, 

& = (k) + w2-  < ,  (*  - 1)  -  »r  <,,(*- 1) 

i  w  -  *L  w + <  (*  ■ - 1)  ■ -  -r 1  <  (a-  - 1) 


and  the  backward  errors  for  any  given  stage  m  are, 

<  (*)  =  (*  “  D +  w?  eL>  W  ~  <  <&_,  (*) 

<  w  =  <*  -  D  -  (*)  +  *7  <*) 


(3-5) 


(3.6) 


(3.7) 
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Figure  11.  T«o  stage  ARMA  lattice  digital  filler. 
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where,  in  Figure  11,  m=  1,2  and  e’a  =  x(k)  and  e}o=y{k).  The  terminal  condition  is 
efm(k)  =  b0  e)m(k)\  m  —  2  .  To  begin  with,  let  b0  equal  unity.  The  initial  conditions  are 
e>0 {k  —  1)  =  0  and  el0{k  —  1)  =  0.  As  with  the  LMS  algorithm,  we  form  an  error  between 
a  desired  signal  d(k)  and  the  output  signal  y(k)  such  that, 

e{k)  =  d{k)-y(k)  (3.8) 

The  instantaneous  gradient  according  to  eq  (3.2)  is  then, 

W)  =  2  e(k)  [  d{k)  -  j )(k)  ]  (3.9) 

Since  the  desired  signal,  d(k),  is  not  a  function  of  the  filter  coefficients,  equation  (3.9) 
reduces  to. 


V(A)  =  2  e(k) 


c\\(k) 


c  -m  ] 


(3.10) 


where  the  quantity 


ov(A) 


C  —  }'(k)  J  is  referred  to  as  the  gradient  estimator.  This  gra¬ 


dient  estimator  must  be  computed  for  each  filter  coefficient  within  the  lattice  structure. 
The  filter  coefficients  are  then  updated  using  the  respective  gradient  estimators.  That  is. 
we  need  to  compute. 


V(k)  «  • 


y(k)  for  /  =  1.2 . 8  and  j=  1,2 . M 


(3.11) 


C  H; 


where  M  is  the  number  of  stages  in  the  ARMA  lattice  filter.  From  equation  (3.5).  the 
gradient  estimates  are  given  by 


cy{k) 


CWi 


ce}x  (k) 
cwi 


+  TT  (*  ~  >)  +  - 


CIV" 


cw. 


d  w J 

-TT  <(*-•) 

CWj 


—  W, 


cebJk  -  1) 


(3.12) 


cvv- 


Let  *k(k)  represent  the  partial  derivatives  of  the  output  y(k)  with  respect  to  the  filter  co¬ 
efficients  and  4>(k)  represent  the  partial  derivatives  of  the  errors  with  respect  to  the  filter 
coefficients.  Using  this  representation  we  can  re-write  equation  (3.12)  as  follows. 
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(3.13) 


'hj  (A)  =  u  (k)  +  6‘i  <  (A-  -  1)  +  wi  ^o,7  (A  —  1)  -  Sli  4o  (A  -  1) 


'  u3  <t>lnU(k  -  1) 


where  <5^,  <5(;,  are  kronecker  deltas  whose  value  is  one  if  and  only  if 
i  —  4  and  j=\  or  /  =  3  and  j  —  1  respectively.  We  compute  ^,,(A)  ,  by  obtaining  re¬ 
cursive  relations  which  calculate  the  partial  derivative  of  the  forward  and  backward  er¬ 
rors  with  respect  to  the  filter  coefficients.  These  relations  are  obtained  from  equations 
(3.6)  and  (3.7)  by  taking  the  partial  derivatives  with  respect  to  the  filter  coefficients, 
namely, 

+U  (*)  =  C,  f j  (*)  +  *"/  (k-\)  +  J?  4>bn_x ,/  (A  -  1)  -  S r/  (A  -  1) 

-<  (A -1) 

<t>l  u  (A)  -  C,  *7  (*)  +  ^+,);  <  (A  -  1)  +  w-r  ’  <  ,J  (A  -  1)  -  <5(3r I)y  <  (A  -  1) 

-*rVZiy(A-l) 


(3.14) 


*7  <*)  "  «7  (*“D  +  *7/  (A)  +  <  c,  17 1*>  -  C  <*> 

*KU  (*)  "  »7  (*"»)-  «?/  (A)  -  C,  U  (A)  +  <5?/  (A) 

+<*Liy<*> 


These  recursive  relations  possess  a  lattice  structure  similar  to  that  of  Figure  11.  with 
delta  components  injected  at  the  summation  nodes.  They  may  also  be  simplified  by  ex¬ 
amining  individual  terms.  Consider  the  general  equation  for  the  forward  error  in  .v.  re¬ 
peated  here  for  continuity. 


£  (*>  =  <*>  +  <  (*  ' ~  1)  -  (A  -  1) 

The  partial  derivative  with  respect  to  each  filter  coefficient  is  expressed  as, 

8efJ  A)  «a_,(A)  cwf  x  „  _  m«?<_i(A-l) 

- - —  — - - - - r  eh  (A  —  1)  +  - ; - 

-  J  ».  /  -.  J  ?n-l  1  '  1  J 

CIV;  CVVj  CtV;  CVV- 

v  ,,  „  ^<_i(a-d 


(3.15) 


(3.16) 


7<_,(A-d-u- 

civ- 
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Since  the  partial  derivative  is  taken  with  respect  to  the  current  filter  coefir’ent  w;  at  time 
A,  the  partial  derivatives  involving  delay  terms  i.e.,  (A  -  1),  are  set  to  ze  This  result 
follows  from  the  realistic  assumption  that  eim  (A:  —  1)  is  a  function  of  u-  (k  —  1)  but  not 
of  vv;  (A)  .  Also  note  that  w;  (A)  is  a  function  of  w;  (A  —  1)  but  not  vice  versa.  With  these 
simplifications  we  reduce  the  equations  of  (3.14)  to, 


K  u  (*>-*/_,  u  «  +  <*?/<.,  <*  - ' )  -  K/  ) 

=  C ,y(*)  +  «+l,/  . '  - 

j*  n.-i  _  i'nj  .X  ( ,s  ,  .. m 


K  u  <*>  -  <"/  (*) + *r  C,  u  <*>  -  ■£< 

k>'  o.i  _  _  j,-x  n,\  i  . 


1) 


n  J  v 


rt  1 7t 


m  ,  v 


and  the  gradient  estimator  is, 

Pij  (*)  =  <PyA  ij  (A)  +  (A  -  1)  -  Sj'  4o  (A  -  1)  (3. 1 S) 


Although  these  are  valid  recursive  relations,  they  are  difficult  to  implement  in  a  lattice 
algorithm.  The  ultimate  goal  is  the  requirement  to  easily  compute  ip,.  (A)  from  the 

available  data.  From  eq  (3. IS)  it  is  evident  that  t//,,  (A)  depends  on  <f>$  (A)  which  in  turn 

requires  knowledge  of  (A)  and  <p}^,,  (A)  but  not  of  (A)  or  (A).  Therefore  the 
three  equations  necessary  to  compute  the  gradient  estimator  are, 

'Pij  W  =  ^  ij  (A)  +  d\{  <  (A  -  1)  -  6\{  c>o  (A  -  1) 

4>lij{k)  =  W  +  *?/<_,  (*  ~  I)  -  *u4m.x  (k  -  0  (3-19) 

4>i  ij  (A)  =  0/M..  i7  (A)  +  ^;+1)y  <  (A  -  1 )  -  d?i+1)JeZm(k-\) 

These  equations  are  dependent  on  the  filter  coefficients  w;\  i~  1.2. 3.4  and 

j  —  1,2 . M  ,  thereby  reducing  by  one-half  the  number  of  computations  required  for 

(A)  and  <j>>/m,,(k).  A  recursive  relation  is  desired  for  p,,(k)  which  does  not  involve 
delta  functions.  Consider  the  four  stage  lattice  filter  with  terminal  condition  b0  equal  to 
unity  such  that  <P)Ml!  (A)  =  (A)  -  The  procedure  for  computing  \p,,  (A)  using  the 

equations  of  (3.19)  is  as  follows: 

1.  Calculate  <p}m0(k)  with  m  equal  to  one  and  letting  i  and  j  range  from  one  to  four. 
Repeat  for  m  —  2,3,4. 

2.  Using  the  terminal  condition  and  expression  for  <p)mtj  (A),  calculate  4>S,VJ  (A). 

This  requires  solving  88  equations,  however,  the  result  is  a  very  simple  recursive  formula 
for  >p,.  (A)  ,  namely. 
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^ij(k)  =  -4  (k-  1)  1=1,3 

r  (3.20) 

'hjW-'ZJk-l)  1=2,4  1 

The  lattice  coefficients  are  calculated  by  substituting  the  recursive  formula  of  eq  (3.20) 
for  the  gradient  estimator  in  equation  (3.3).  The  adaptive  coefficient  update  equation  is, 

(*)  =  w{  (k-l)-fi  e(k)  4>ij  (k)  (3.21) 

Since  the  gradient  estimator,  and  therefore  the  gradient,  is  the  same  for  vv>;  ,  /=  1,3  and 
w;,  i  —  2,4  ,  it  follows  that  wj  =  w;  and  iw)  =  wi.  The  number  of  filter  coefficients  re¬ 
quired  to  update  the  lattice  filter  is  reduced  to  two,  i.e.,  vvj  and  wj.  Furthermore,  from 

the  symmetry  of  the  lattice  structure,  the  following  equalities  between  filter  coefficients 
are  assumed. 


=  wy 


xv3 


”4 


u> = 


w8 


(3-22) 


Incorporating  these  equalities  with  those  derived  by  the  gradient  estimator  produces  the 
elementary  ARMA  lattice  section  of  Figure  12  where. 

U-2  =  u-j  =  \vJ5 

j  j  j 

u/j  =  w'  =  u/7 


=  ui 


=  vvi  =  k, 


(3.23) 


To  prove  that  these  coefficient  reductions  are  valid,  a  computer  generated  solution  using 
the  Fortran  program  of  Appendix  D  was  compared  to  hand  analysis  of  a  second  order 
transfer  function  and  lattice  filter.  The  output  of  the  ARMA  digital  lattice  filter  was  first 
put  into  difference  equation  form  and  then  compared  to  the  known  transfer  function. 
From  this  comparison,  lattice  coefficients  were  computed.  Details  of  this  analysis  are  as 
follows.  Consider  a  two  stage  ARMA  digital  lattice  filter  comprised  of  the  reduced  ele¬ 
mentary  section  shown  in  Figure  13  and  a  transfer  function  of  the  form. 


bo  +  b\  z  bj  z 

1  +  fl]  Z  +  tJj  Z 


(3.24) 
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Figure  12.  Simplified  elementary  ARMA  lattice  section. 

The  output  of  the  lattice  filter  can  be  written  in  difference  equation  form  by  carrying  out 
the  following  steps:  (i)  start  with  the  output  of  the  lattice  filter  equation  (3.5),  (ii)  sub¬ 
stitute  expressions  for  the  forward  and  backward  errors,  equations  (3.6)  and  (3.7).  re¬ 
spectively.  into  equation  (3.5)  and  (iii)  carry  out  the  algebra.  A  detailed  derivation  of  this 
difference  equation  is  given  in  Appendix  E.  The  difference  equation  in  its  final  form  is 
given  by, 

y(k)  =  x(k)  +  2(  ivj  +  w  '  +  w2  w2  )x(k  —  1)  +  2  w2  x{k  —  2) 

— 2(  wj  +  w2  w\  -t-  VV|  wf  M*  "  1)  -  2  w?.v(*  -  2) 

which  can  be  written  in  the  transfer  function  form  as, 


Hiz) 


12.  1  2  x  —  1  .  2  —2 

1  +  2(  IV2  +  VV,  IV,  +  VV,  VV2  )  2  +  2  W,  2 

.  .  i  !  i  2  !  i  2  x  -l  "  I  2 

1  +  2(  VV,  +  u'2  Wj  +  VV,  VV,  )  2  +2  VC,  2 


(3.26) 


Comparing  the  lattice  filter  transfer  function,  H(z)  with  the  known  filter  transfer  func¬ 
tion  Hj(z),  produces  the  following  relationships  between  filter  coefficients, 


r  1  |  12.  1  2  \ 

bx  =  2(  w2  +  wx  wx  +  vv2  vv2  ) 

a{  =  2(  vv,1  +  vvj  w|  4-  vv,  wf  ) 

b2  —  2  vv2 

,  2 
a2  =  2  vv. 


(3.27) 


Solving  for  vv^  and  vvf  in  terms  of  the  known  transfer  function  coefficients  b2  and  a2  and 
then  substituting  these  results  into  the  expressions  for  4>,  and  au  respectively,  yields  the 
following. 


bx  =  a2  wj  +  (2  +  b2)  vv] 
tf,  =  (  2  4-  a2  )  vv,1  +  b2  w2 


or  in  matrix  form, 


a2  2  4-  b2 
2  +  a2  b2 


vv. 


vv, 


(3.2S) 


(3.29) 


and  solving  for  the  lattice  coefficients,  we  have 


*  1 

M'l 

1 

b2 

-  (2  +  b2)  ’ 

b\ 

1 

. u'2  - 

a2  h2  -  (2  4  b2)  (2  4-  a2) 

_  —  (2  -f  a2) 

a2 

.  . 

and 


(3.31) 


Now  that  a  method  of  converting  between  lattice  and  transfer  function  coefficients  for 
a  second  order  system  has  been  established,  we  consider  the  specific  transfer  function 
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H/z) 


1  -  O.S  r~‘  +  1.78  z~2 
1  -  0.S9  z-1  +  0.25  z~2 


(?-32) 


where 


b0  =  1.0  bx  =  —0.80  A2=1.78 

=  -0.89  a2  =  0-25 

From  equations  (3.30)  and  (3.31)  the  lattice  coefficients  are  calculated  as, 

wf  =  0.125,  w'2  =  0.890,  w}  =  -0.240719,  w2'  =  -0.195719 

Values  for  the  steady-state  lattice  coefficients  were  computed  using  the  Fortran  program 
in  Appendix  D  and  are  shown  below.  Convergence  aspects  of  both  the  lattice  coeffi¬ 
cients  and  output  error  are  shown  in  Figure  13. 

w?  =0.124982,  wl  =0.890003,  =  -0.240710,  wj  =-0.195711 


these  results  confirm  the  validity  of  the  derived  adaptive  lattice  algorithm  and  the  design 
of  a  new  elementary  lattice  section  shown  in  Figure  12. 

The  current  adaptive  lattice  algorithm  assumes  that  the  terminal  condition  is  unity. 
This  is  generally  not  the  case  in  practice.  We  now  extend  this  adaptive  algorithm  to  the 
more  general  case  where  the  terminal  condition  is  an  arbitrary  constant.  The  recursive 
relation  which  updates  b0  is  similar  to  those  which  update  the  other  lattice  filter  coeffi¬ 
cients.  The  update  equation  for  *0  is  given  by 


b0{k)  =  b0(k-  1 )~ne{k) 


ce(k ) 
cb0(k) 


(3.33) 


ce(k) 


The  gradient  estimator  ■  ■  ,  is  calculated  using  equations  (3.5),  (3.8)  and  the  fact  that 
the  desired  signal  d'kj  is  not  dependent  on  b0.  The  gradient  estimator  for  b0  is  written 
as. 


ii-m  *#*) 


bbn 


Cbn 


■  +  ' 


1 

cwA 

cbn 


x  ,  <?<(*-!) 
•»*-'> + -4  “V— 


-  1 
Cbn 


—  w. 


<(*- 1) 

~  0 

Cbn 


(3.34) 


Since  the  partial  derivative  is  taken  with  respect  to  *0  at  time  k,  this  reduces  to, 
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cy(k)  ce}(k) 

cb0  cb0 


Similarly, 

i  <*>  -  w + -r 1  <(*-D-  «r 1 1  «*  - » 

and  the  partial  derivative  with  respect  to  b0  is 

K  (*)  _  *L,  (*) 

c60  Cib0 


The  terminal  condition  is, 


&(*)«**>  <£,(*) 


(3.36) 


(3.37) 


(3.3S) 


Taking  the  partial  derivative  of  equation  (3.3S)  with  respect  to  b0  yields  (k),  and  the 
recursive  equation  to  update  ba  ( k )  becomes, 

b0(k)  =  b0(k-l)-ne(k)efM(k)  (3.39) 


The  gradient  estimators  for  the  lattice  filter  coefTicients  are  scaled  by  the  arbitrary  con¬ 
stant  b0,  since  at  the  terminal  condition  {k)  =  60  4>}m,j  {k)  and  b0  is  propagated 
through  the  calculations.  The  gradient  estimators  become, 

'I'tjW  =  ~b0  el  (k  -  1 )  /-lt3 

x  1  (3.4d) 

i>ij  (k)  =  b0  c£_(  (k  —  \)  i  =  2,4 


To  test  this  more  general  adaptive  lattice  algorithm  the  output  of  a  known  transfer 
function  with  ba  equal  to  0.5  was  compared  to  the  output  of  the  ARMA  digital  lattice 
filter.  The  second  order  transfer  function  used  was, 


H/z)  - 


0.5  —  0.4  z~‘  4-  0.89  z~2 
1  —  0.89  z-1  +  0.25  z~2 


(3.41) 


The  computer  generated  steady  state  lattice  coefficients  are  given  below  and  convergence 
aspects  shown  in  Figure  14. 

b0  =  0.499946,  =  -0.120428.  w2  =  0.593237,  »•’  -  -0.447423,  uf  -  0. 166320 
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In  order  to  maintain  the  same  adaptive  time  constant  and  misadjustment  at  each 
stage  in  the  lattice,  the  convergence  constant  is  normalized  by  the  power  level  at  each 
stage  [Ref,  15].  Therefore,  we  can  write  equations  (3.21)  and  (3.39)  as, 


J,  (k)  =  w/(k-  1) - r— 

Oj  (A) 


b0(k)  =  b0(k-  1) - f— 

y  (A) 


*(A)  tij  (A) 

e(k)  efst  (A) 


(3.42) 


where  n  is  the  convergence  constant  and  a)  ( k )  and  yJ  (A)  are  estimates  of  the  power  at 
the  j,h  stage  for  w;  and  b0  respectively  and  computed  as  follows: 


Oj  (A)  =  p  oj  (k  —  1)  +  (1  —  p)  tfj(k) 

>’2  (A)  =  P  y2  (A  —  1)  +  (1  —  p)  [  (A)  ]2 


(3.43) 


Writing  equation  (3.42)  using  the  notation  adopted  for  the  reduced  elementary  ARMA 
lattice  section  we  obtain 


rj  (A)  =  rj  (k  -  1) - f—  e(k)  4  (A  -  1) 

Oj  (A) 

kj  [k)  =  kj  (k  -  1)  -  -yJ—  e(k)  4  (k  -  1) 
Oj  (A) 

MA)»MA  — 1) - r~  e(k)efx  (k) 

7  (A) 


(3.44) 


In  the  above  equations  p  is  a  weighting  parameter,  0  <,  p  <  1,  which  distributes  the 
amount  of  weight  given  the  past  power  level  or  current  sample.  Normalized  convergence 
constants  are  used  in  all  examples  of  this  thesis.  The  adaptive  lattice  algorithm  is  sum¬ 
marized  in  Table  1. 

In  summary  ,  we  have  derived  an  adaptive  algorithm  based  on  the  LMS  theory  of 
adaptive  coefficient  computation.  This  new'  adaptive  algorithm  easily  updates  the  lattice 
coefficients  by  using  available  data.  The  original  requirement  to  update  eight  coefficients 
of  an  elementary  ARMA  lattice  section  was  reduced  to  updating  only  two  coefficients 
and  still  being  able  to  describe  the  lattice.  The  algorithm  is  general  in  that  it  applies  to 
systems  whose  terminal  condition  is  an  arbitrary  constant.  The  validity  of  this  algorithm 
was  demonstrated  through  comparisons  between  hand  analysis  and  computer  simu¬ 
lation.  In  the  next  chapter,  we  further  demonstrate  the  convergence  of  this  algorithm. 
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Table  1.  SUMMARY  OF  ADAPTIVE  LATTICE  ALGORITHM 


With  given  initial  conditions,  ef0  ~  x{k),  e-}g  =  >(A),  e>B  {k  -  1)  =  (A  -  1)  =  0  and  all 
lattice  coefficients  zero, 

Step  1:  for  k=  l,m  compute 

<£(*)  -  eL,  <*>  +  “i"  (*-!)- (*  -  1) 

£  <*>  -  «L,  (*)  +  <(*-')-  "C  <,  <*  - 1) 

with  output  e*  (A) 

Step  2:  for  k=  l,m  compute 

<  <*>  -  (*-d+ <£.,  (*)  -  <  <L,  w 

<  <*)  -  <*-*>-  <*> + <£_,  <*> 

Step  3:  Update  coefficients 


r,(A)  =  r;(A-l) - e(A)  (A  -  1) 

(A) 

A,  (A)  =  kj(k  -  1)  —  c(A)  ^  (A  -  1) 
°j  (A) 

W  =  ^-D — rrr^)e/  (A) 
y  (*) 


Step  4:  Repeat  for  next  iteration  i.e.  return  to  step  1. 


IV.  EXPERIMENTAL  RESULTS 

The  adaptive  lattice  algorithm  derived  in  Chapter  III  is  now  computer  simulated  to 
study  its  convergence  performance.  The  system  identification  mode  of  adaptive  filtering 
is  considered  for  this  purpose.  Figure  15  shows  a  general  system  identification  config¬ 
uration.  The  systems  considered  are  time-invariant  and  linear.  Notice  that  we  apply  the 
same  input,  white  noise  in  general,  to  both  the  reference  system  and  the  adaptive  lattice 
filter  which  is  modeling  the  system.  The  criterion  in  this  configuration  is  to  minimize  the 
mean-squared  error  between  the  system  and  filter  outputs.  Thus,  in  this  context,  the 
adaptive  algorithm  continuously  updates  the  lattice  filter  parameters  in  order  to  mini¬ 
mize  the  mean-squared  error. 

The  adaptive  algorithm  is  realized  as  summarized  in  Table  1.  As  we  mentioned  in 
Chapter  III.  the  two  important  parameters  of  the  algorithm  are  the  adaptation  constant 
H  and  the  weighting  constant  p  .  In  what  follows,  we  shall  consider  convergence  studies 
of  both  second  and  third  order  reference  systems  (fixed  filter  transfer  functions).  Con¬ 
sider  the  following  reference  system  with  transfer  function, 

1  4-  0.2  z~‘  -  0-35  z~2 
1  -  1.4  z-1  4-  0.85  z~2 

« 

This  system  has  complex  poles  and  simple  zeros  located  at  z  =  (0.7  ±j0. 6)  and 
z  =  0.5, —0.7,  respectively.  Using  a  convergence  constant  p  =  0.01  and  power  level 
weighting  factor  p  =  0.45,  the  adaptive  ARMA  digital  lattice  filter  which  models  the 
above  system  has  the  following  steady-state  lattice  parameters, 

terminal  condition  b0  =  0.999202 

lattice  coefficients  r,1  =  0.352580 
r}  =  -0.174355 
k}  =  -0.448228 
kf  =  0.425060 

Convergence  properties  of  the  lattice  coefficients  and  error  are  shown  in  Figure  16.  The 
mean-squared  error  was  minimized  after  approximately  1700  iterations  at  which  time  the 
lattice  coefficients  reached  their  steady-state  values.  When  the  value  of  the  convergence 
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constant  p  was  modestly  increased,  convergence  degraded  rapidly.  Also,  when  the 
weighting  factor  p  was  increased,  convergence  deteriorated  quickly.  From  this,  we  con¬ 
clude  that  the  convergence  constant  is  the  more  sensitive  input  parameter. 

Let  us  consider  another  second  order  dynamic  system  with  transfer  function, 


0.5-  0.2  -  0.4-15  :~2 

1  -z_i  +0.9-1  z-2 


This  system  has  complex  poles  at  z  =  (0.5  ±jO.S)  and  complex  zeros  located  at 
z  =  (0.2  ±y0.7).  Using  a  convergence  constant  p  =  0.005  and  power  level  weighting 
factor  p  —  0.97,  the  adaptive  ARMA  digital  lattice  filter  which  models  this  system  has 
steady-state  lattice  parameters  as  follows, 


terminal  condition  60  =  0.500027 


lattice  coefficients  r,’  =  0.104654 
r]  =  0.296658 
k\  =  -0.429232 
k]  =  0.627718 
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Convergence  properties  of  this  adaptive  filter  are  shown  in  Figure  17.  In  this  example, 
convergence  was  obtained  after  approximately  2500  iterations.  Again,  by  changing  the 
values  of  p  and  p  slightly,  covergence  deteriorated  with  the  convergence  constant  p  being 
the  more  sensitive  parameter. 

Next  we  consider  a  third  order  reference  system  with  known  transfer  function, 


///;)  = 


0.5  -  0.95  2~'  4-  1 .33  z~2  -  0.979  z~3 
1-1.69  z_1  +  0.962  z-2  -  0.2  z-3 


The  adaptive  ARMA  digital  lattice  filter  which  describes  this  system  has  the  following 
steady-state  lattice  parameters, 

terminal  condition  b0  =  0.499970 

lattice  coefficients  r!  =  —0.328447 
r]  =  0.399472 
r\=  -0.652706 
k\  =  — 0.S2173S 
k]  =  -0.091111 
k\  =  -0.133333 


These  parameters  were  obtained  using  a  convergence  constant  p  =  0.015  and  power  level 
weighting  factor  p  =  0.9.  Convergence  properties  of  this  adaptive  filter  are  shown  in 
Figure  IS.  Steady-state  values  for  the  lattice  coefficients  were  obtained  after  approxi¬ 
mately  7100  iterations.  It  is  reasonable  to  assume  that  a  third  order  system  will  converge 
more  slowly  than  a  second  order  system.  The  number  of  iterations  required  for  this  third 
order  system  to  converge  is  consistant  with  convergence  rates  of  other  adaptive  algo¬ 
rithms  which  model  third  order  systems  [Ref.  16].  The  input  parameter  p  was  again 
found  to  be  the  more  sensitive  parameter. 

In  all  the  previous  examples,  the  values  of  p  and  p  may  or  may  not  be  optimum 
values.  That  is,  an  exhaustive  search  of  all  combinations  of  p  and  p  was  not  performed 
to  demonstrate  convergence  of  the  algorithm.  Nevertheless,  a  number  of  different  ways 
of  realizing  the  value  of  the  convergence  constant  p  have  been  reported  in  the  literature. 
In  one  method,  Mikhael  et.  al.  [Ref.  17]  have  obtained  a  variable  p  by  using  a  self  opti¬ 
mizing  technique.  In  this  method,  p  is  calculated  from  the  input  data  as  an  iteration 
process  and  is  individually  determined  for  each  filter  parameter.  In  another  method  p  is 
chosen  by  using  a  variable  step  LMS  technique  [Ref.  IS],  where  the  range  of  p  is 
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ERROR  MAGNITUDE 
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Figure  17. 


Second  order  ARMA  lattice  Otter,  terminal  condition  of  0.5. 
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specified  by  pmix  and  ^min  which  are  within  the  bounds  described  by  equation  (3.4)  of 
Chapter  III.  These  techniques  of  choosing  p  during  the  adaptation  process  have  been 
shown  to  improve  filter  convergence.  They,  however,  require  additional  computations 
to  achieve  this  faster  convergence.  When  a  combination  of  the  parameters  p  and  />  was 
obtained  which  yielded  convergence,  these  values  were  chosen  for  examples.  Besides  the 
examples  reported,  simulation  studies  have  been  carried  out  for  several  other  cases.  In 
all  cases,  however,  definite  convergence  of  the  algorithm  has  been  observed. 

In  summary,  we  have  demonstrated  through  computer  simulation  that  the  derived 
adaptive  lattice  algorithm  is  suited  for  system  identification  modeling.  Furthermore,  we 
have  shown  that  there  is  flexibility  in  choosing  the  values  of  the  convergence  constant 
p  and  weighting  factor  p  .  Some  techniques  for  selecting  (computing)  the  value  of  the 
convergence  constant  have  been  introduced.  These  methods  improve  convergence  at  the 
cost  of  additional  computations. 

A.  CONCLUSIONS 

In  this  thesis  we  have  demonstrated  that  the  ARMA  parameter  estimation  algo¬ 
rithm  proposed  in  [Ref.  4:  pp.  619-621]  is  a  valid  method  for  obtaining  approximations 
to  reference  models.  Furthermore,  the  criterion  used  to  derive  the  algorithm  is  a  gener¬ 
alized  form  of  the  Mullis-Roberts  criterion  for  least  squares  modeling.  The  AR  and  MA 
parameters  of  the  ARMA  model  can  be  updated  independently  as  their  respective  orders 
increase  by  one.  From  the  recursive  prediction  error  formulas,  an  ARMA  digital  lattice 
filter  was  designed  with  arbitray  AR  and  MA  orders. 

For  the  ARMA  digital  lattice  filter,  we  derived  an  adaptive  lattice  algorithm.  This 
algorithm  was  based  on  the  least  mean  square  method  of  optimizing  coefficients.  The 
derived  adaptive  lattice  algorithm  can  easily  compute  the  values  of  the  lattice  coefficients 
from  available  data.  The  algorithm  simplified  the  number  of  coefficents  required  to  be 
updated  from  eight  coefficients  per  elementary  lattice  section  to  only  two  such  that  the 
filter  can  be  completely  described.  This  savings  in  computational  effort  makes  the  al¬ 
gorithm  attractive  for  identification  of  unknown  systems  since  many  systems  require  an 
ARMA  model  for  parsimonious  modeling. 

Convergence  of  the  adaptive  lattice  algorithm  was  demonstrated  with  several  exam¬ 
ples  in  Chapter  III.  The  number  of  iterations  required  before  convergence  varied  greatly 
between  second  and  third  order  models  as  well  as  within  second  order  models.  Optimum 
convergence  rates  were  not  sought  after  as  much  as  proving  the  convergence  of  the  al¬ 
gorithm.  Rapid  convergence  rates  were  demonstrated  in  Chapter  II  for  a  second  order 
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system  upon  completion  of  an  extensive  search  for  the  optimun  values  of  the  conver¬ 
gence  constant  and  power  level  weighting  factor. 

Although  a  method  of  converting  between  direct  form  and  lattice  realizations  for 
second  order  ARM  A  filters  was  developed,  a  general  algorithm  to  perform  this  trans¬ 
formation  given  the  ARMA  lattice  filter  design  was  not  obtained.  When  solving  for  a 
transformation  between  filter  realizations  of  third  order  the  solution  is  hindered  by 
nonlinearities. 

The  objectives  of  the  thesis  were  sucsessfully  accomplished.  Some  suggestions  for 
future  work  include  the  following:  (i)  extensive  theoretical  analysis  for  determining  op¬ 
timum  values  for  n  ,  (ii)  derivation  of  a  generalized  algorithm  which  converts  any  given 
ARMA  transfer  function  into  a  set  of  lattice  parameters,  (iii)  development  of  theoretical 
convergence  models  for  the  ARMA  adaptive  lattice  algorithm  and  analysis  of  these 
models  and  (iv)  application  of  the  adaptive  lattice  filter,  both  analysis  and  synthesis 
forms,  in  modeling  such  practical  signals  as  speech.  ARMA  lattice  filter  modeling  has 
considerable  application  potential  because  of  its  very  accurate  modeling  of  nearly  any 
signal  or  system  of  interest. 
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APPENDIX  A.  MAIN  PROGRAM  TO  ESTIMATE  ARMA  PARAMETERS 


C  THIS  PROGRAM  COMPUTES  THE  ARMA  PARAMETERS  AS  THE  AR  OR  MA  ORDER 

C  OF  AN  ARMA  MODEL  INCREASES  BY  ONE.  IT  USES  THE  ARMA  PARAMETER 

C  ESTIMATION  ALGORITHM  PROPOSED  BY  MIYANAGA,  NAGAI  AND  MIKI. 

C 

C  VARIABLE  DEFINITIONS 

C 

C  VN  -  INPUT  VECTOR  CONTAINING  DATA  GENERATED  AS  WHITE  GAUSSIAN 

C  NOISE 

C  Y  OUPUT  VECTOR  OF  DIFFERENCE  EQUATION  WITH  VN  AS  INPUT 

C  RX  -  AUTOCORRELATION  DATA  OF  INPUT  VN 

C  RY  AUTOCORRELATION  DATA  OF  OUTPUT  Y 

C  RXY  -  CROSSCORRELATION  DATA  OF  INPUT  AND  OUTPUT 

C  RYX  -  CROSSCORRELATION  DATA  OF  OUTPUT  AND  INPUT 

C  NDATA  -  NUMBER  OF  INPUT  DATA  POINTS 

C  KDATA  -  NUMBER  OF  INITIAL  DATA  POINTS  TO  DISREGARD 

C  XI  -  INPUT  DATA  VECTOR  AFTER  DISCARDING  KDATA  POINTS 

C  Y1  OUTPUT  DATA  VECTOR  AFTER  DISCARDING  KDATA  POINTS 

C  A  -  VECTOR  CONTAINING  CALCULATED  AR  PARAMETERS 

C  B  VECTOR  COINTAINING  CALCULATED  MA  PARAMETERS 

C  TXA  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 

C  INPUT 

C  TYA  -  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 
C  OUPUT 

C  TXB  -  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 
C  INPUT 

C  TYB  -  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 
C  OUTPUT 

C  GA  VECTOR  CONTAINING  COEFFICIENTS  FOR  BACKWARD  PREDICTION  OF 

C  INPUT 

C  GB  VECTOR  CONTAINING  COEFFICIENTS  FOR  BACKWARD  PREDICTION  OF 

C  INPUT 

C  ZTA  -  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 
C  OF  OUTPUT  VALUES 

C  ZTB  -  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 
C  OF  OUTPUT  VALUES. 

C  NI  LENGTH  OF  DATA  VECTOR  XI 

C  NK  -  LENGTH  OF  DATA  VECTOR  XI  MINUS  ONE  USED  TO  START 

C  CORRELATION  COMPUTATIONS. 

C  NS  DESIRED  AR  ORDER  OF  ARMA  MODEL. 

C  NT  DESIRED  MA  ORDER  OF  ARMA  MODEL 

C  KS  CURRENT  AR  ORDER  OF  UPDATE 

C  KT  -  CURRENT  MA  ORDER  OF  UPDATE 

C  VX  EXPECTED  VALUE  OF  PREDICTION  ERROR  FOR  INPUT  SQUARED. 

C  VY  EXPECTED  VALUE  OF  PREDICTION  ERROR  FOR  OUTPUT  SQUARED. 

C  VXY  EXPECTED  VALUE  OF  PRODUCT  OF  PREDICTION  ERROR  FOR  INPUT 


n  o  o  o  o  Kj  i-*  ooo 


C 

C 

c 

c 

c 

c 

c 


AND  OUTPUT. 

VG  -  EXPECTED  VALUE  OF  BACKWARD  PREDICTION  ERROR  OF  INPUT 
SQUARED. 

VZ  -  EXPECTED  VALUE  OF  BACKWARD  PREDICTION  ERROR  OF  OUTPUT 
SQUARED. 

VGZ  -  EXPECTED  VALUE  OF  PRODUCT  BETWEEN  BACKWARD  PREDICTION 
ERRORS  OF  INPUT  AND  OUTPUT. 

DIMENSION  VN( -5: 2006) ,Y(-5: 2006) ,RX(0: 2000) ,RY(0:  2000) ,RXY(0:  2000) 
DIMENSION  RYX(0: 2000) ,X1(2000) ,Y1(2000) ,X( 12) ,A(0: 21) ,B(0: 21) 
DIMENSION  TXA(0:  21) ,TYA(0:  21) ,TXB(0:  21) ,TYB(0:  21) ,GA(0:  21) 
DIMENSION  GB(0: 21),ZTA(0: 21) ,ZTB(0: 21) 

INPUT  DATA  INFORMATION 

WRITE  (6,1) 

FORMAT  (/'  ENTER  THE  NUMBER  OF  DATA  POINTS:  ’) 

READ  (6,*)  NDATA 
WRITE  (6,2) 

FORMAT  (/'  ENTER  NUMBER  OF  INITIAL  DATA  POINTS  TO  DISREGARD:  ') 

READ  (6,*)  KDATA 

*********^V*******^V***^********Vr***********>V****'>'r**iVTV*Vr********** 


INITIALIZE  ARRAYS 


DO  10  L=-5, NDATA 
VN(L)=0 
Y(L)  =0 
10  CONTINUE 

DO  15  L=0 , 20 
A(L)=0 
B(L)=0 
TXA(L)=0 
TXB( L)=0 
TYA(L)=0 
TYB( L)=0 
GA(L)  =0 
GB(L)  =0 
ZTA(L)=0 
ZTB(L)=0 
15  CONTINUE 

DO  20  L=0, 1999 
RX(L)=0 
RY(L)=0 
RXY(L)=0 
RYX(L)=0 
20  CONTINUE 

DO  25  L=1 , 12 
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X(L)=0 

25  CONTINUE 

C 

C 

C  GENERATE  WHITE  GAUSSIAN  NOISE  INPUT 

C 

ISIZE=NDATA 

IX=152255 

ISORT=0 

MUL=2 

DO  30  K=1 , ISIZE 

CALL  SRND( IX,X, 12,MUL, ISORT) 

XT=-6. 0 

DO  35  1=1,12 
35  XT=XT+X( I) 

VN(K)=XT 
30  CONTINUE 

C 

C 

C  COMPUTE  OUTPUT  OF  REFERENCE  MODEL  FILTER  AND  DISREGARD  SPECIFIED 
C  NUMBER  OF  DATA  POINTS 

C 

DO  40  L=1 ,NDATA 

C  Y(L)=VN(L)+1.  6*Y(L-l)-0.  95*Y(L-2) 

C  Y(L)=VN( L)+0.  2*Y(L-l)-0.  62*Y(L-2)+0. 152*Y(L-3)-0.  3016*Y(L-4) 

C  Y(L)=VN(L)-0.  2*VN(L-l)+0.  62*VN(L-2)-0.  152*VN(L-3)+0.  3016*VN(L-4) 

C  Y(L)=VN(L)+0.  2*VN(L-1) -0.  99*VN(L-2)+0. 2*Y(L-l)-0.  62*Y(L-2)+0.  152*Y 

C  &(L-3)-0. 3016*Y(L-4) 

C  Y(L)=VN(L)-1.  6*VN(L-1)+1. 45*VN(L-2)+l. 2*Y(L-l)-0.  72*Y(L-2) 

C  Y(L)=VN(L)+0.  2*VN(L-l)-0.  35*VN(L-2)  +  l.  4*Y(L-l)-0.  85*Y(L-2) 
Y(L)=0.5*VN(L)-0.  2*VN(L-l)+0.  445*VN(L-2)+Y(L-l)-0.  94*Y(L-2) 

C  Y(L)=VN(L)-2.  7*VN(L-l)+3. 21*VN(L-2) -1. 595*VN(L-3)+l.  95*Y(L-1)-1.  62 
C  &*Y(L-2)+0.  54*Y(L-3) 

C  Y( L)=VN( L) -1. 0*VN(L-l)+0. 89*VN(L-2)+0.  40*Y(L-l)-0.  2121*Y(L-2) -0.  20 

C  &894*Y(L-3)-l.  810373*Y(L-4) 

C  Y(L)=0.  5*VN(L) -0.  95*VN(L-1)+1. 33*VN(L-2)-0. 979*VN(L-3)+l. 69*Y(L-1) 

C  &-0. 962*Y(L-2)+0.  20*Y(L-3) 

C  Y(L)=0.  5*VN(L)-0.  4*VN(L-l)+0.  89*VN(L-2)+l.  69*Y(L-l)-0.  962*Y(L-2)+0 

C  &. 2*Y(L-3) 

C  Y(L)=0.  5*VN(L)-0.4*VN(L-1)+0. 89*VN(L-2)+0.  2*Y(L-l)+0.  25*Y(L-2)-0.  0 

C  &5*Y(L-3) 

C  Y(L)=0.  5*VN(L) *0.  4*VN(L-l)+0. 89*VN(L-2)+0. 89*Y(L-l)-0.  25*Y(L-2) 

C  Y(L)=.  0154*VN(L)+.  0642*VN(L-l)+0.  0642*VN(L-2)+0. 0154*VN( L-3)+l.  99* 

C  &Y(L-1) -1.  57*Y(L-2)+0.4583*Y(L-3) 

C  Y(L)=0.  5*VN(L)+0.  256*VN(L-l)+0.  1234*VN(L-2)+0.  0987*VN(L-3) 

40  CONTINUE 
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LJ=KDATA+1 
DO  45  L=LJ , NDATA 
LK=L-KDATA 
Xl( LK)=VN( L) 

Y1(LK)=Y(L) 

45  CONTINUE 

C  WRITE  (*,77)  (Y1(K) ,  K=200,1800) 

77  FORMAT  (5( 1X.F10. 6)) 

C 

C 

C  COMPUTE  AUTO-CORRELATION  AND  CROSS -CORRELATION  TERMS 

C 

46  NI=NDATA-KDATA 
NK=NI-1 

CALL  CORREL  (NI ,50 ,X1 ,Y1 ,RX,RY ,RXY,RYX,NK) 

47  DO  50  L=0 , 10 

WRITE  (*,200)  RX(L) ,RY(L) ,RXY(L) ,RYX(L) 

WRITE  (9,200)  RX(L) ,RY(L) ,RXY(L) ,RYX(L) 

WRITE  (9,201) 

50  CONTINUE 

WRITE  (9,201) 

200  FORMAT  (2X,4(2X,F14. 9)) 

201  FORMAT  (’  ’) 

C 

Q  **Vf**Vr****,V*********^ViV******iV*Vf*i^*TV************iWr*TV*ilr>V**'iV***->V->V** 

c 

c  INPUT  THE  DESIRED  AR  AND  MA  ORDERS  THEN  DEFINE  INITIAL  CONDITION’S 

C 

48  WRITE  (6,3) 

3  FORMAT  (/'  ENTER  THE  DESIRED  AR  ORDER: ’) 

READ  (6,*)  NS 

WRITE  (6,4) 

4  FORMAT  (/'  ENTER  THE  DESIRED  MA  ORDER: ') 

READ  (6,*)  NT 

C 

KS=0 

KT=0 

A(0)=1. 0 
VX=RX( 0 ) 

VY=RY(0) 

VXY=-RYX(0) 

VG=RX( 0) 

VZ=RY( 0) 

VGZ=-RYX(0) 

TXB( 0)  =  1.  0 
TYA(0)=1.  0 
GB(  0)  =1.0 
ZTA(0)=1.  0 
C 
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c 

C  ESTIMATE  THE  ARMA  PARAMETERS 
C 

300  IF  (NT.  EQ.O.  AND.KS.LT.  NS)  THEN 

KS=KS+1 

CALL  NEWAR( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , 2TB 
& , RX , RY , RXY , RYX ,  A ,  B ) 

GOTO  300 
ELSE 

301  IF  (NS. EQ. 0.  AND.  KT.  LT.  NT)  THEN 

KT=KT+1 

CALL  NEWMA( KS ,  KT ,  VX ,  VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
&,RX,RY,RXY,RYX,A,B) 

GOTO  301 
END  IF 
ENDIF 
C 

IF  (NS.NE.  O.OR.  NT.  NE.  0)  THEN 

302  IF  (NS.  GE.  NT.  AND.  NT.  NE.  0.  AND.  KT.  LT. NT)  THEN 

KS=KS+1 

CALL  NEWAR( KS ,  KT ,  VX ,  VY , VXY , VG , VZ , VGZ , TXA , TXB ,  TYA , TYB , GA , GB , ZTA , ZTB 
&,RX,RY,RXY,RYX,A,B) 

KT=KT+1 

CALL  NEWMA( KS ,  KT ,  VX ,  VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
& , RX , RY , RXY , RYX , A , B ) 

GOTO  302 
ELSE 

303  IF  (KS. LT. NS)  THEN 

KS=KS+1 

CALL  NEWAR(KS ,KT, VX, VY, VXY,VG,VZ, VGZ,TXA,TXB ,TYA,TYB ,GA,GB , ZTA, ZTB 
& , RX , RY , RXY , RYX , A , B ) 

GOTO  303 
ENDIF 
ENDIF 
ENDIF 

n 

o 

C  **VrVc*Vc**Vf*Vr*VVyf**VfVc**TV******Vr*****>Vyr*>VVf^Vryr****ynVVf*Vr'jV******’5V*Vf**Vt-yryrVr 

c 

C  PRINT  ESTIMATED  ARMA  PARAMETERS 

C 

WRITE  (*,211) 

WRITE  (9,211) 

WRITE  (*,210)  (A(K) ,  K=1,KS) 

WRITE  (9,210)  ( A(K) ,  K=1,KS) 

WRITE  (*,211) 

WRITE  (9,211) 

211  FORMATC  ') 

WRITE  (*,210)  (B(K) ,  K=0,KT) 
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o  o  o  u  o 


WRITE  (9,210)  (B(K) ,  K=0,KT) 
210  FORMAT  ('  ’ , 1X,4(2X,F13. 10)) 
STOP 
END 


SUBROUTINE  TO  COMPUTE  CORRELATION  TERMS 

SUBROUTINE  CORREL( N , LAG , X , Y , RX , RY , RXY , RYX , NK1 ) 

REAL  X(0: NK1) ,  Y(0:  NK1) ,RX(0:  2000), RY(0:  20 00),RXY(0:  2000) 
REAL  RYX( 0:  2000) , SUM1 , SUM2 , SUM3 , SUM4 
DO  70  K=0 , LAG 
NJ=N-1-K 
SUM 1=0 
SUM2=0 
SUM3=0 
SUM4=0 
ANK=NJ 

DO  60  J=0 ,NJ 

SUM1=SUM1+X( J+K)*X( J) 

SUM2=SUM2+Y(J+K)*Y(J) 

SUM3=SUM3+X( J+K )*Y( J ) 

SUM4=SUM4+X( J)*Y( J+K) 

60  CONTINUE 

RX(K)=SUM1/ANK 
RY ( K ) =SUM2 / ANK 
RYX (K)=SUM3/ ANK 
RXY(K)=SUM4/ANK 
70  CONTINUE 
RETURN 
END 
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APPENDIX  B.  SUBROUTINE  FOR  MAIN  PROGRAM 


SUBROUTINE  NEWAR(KS,KT,VX,VY,VXY>VG,VZ>VGZJTXA)TXB,TYA,TYB,GA,GB,Z 
&TA , ZTB , RX , RY , RXY , RYX , A , B ) 

C 

C  THIS  SUBROUTINE  COMPUTES  AR  PARAMETER  VALUES  FOR  AN  ARMA  MODEL  AS 
C  THE  AR  ORDER  INCREASES  BY  ONE. 

C 

C  VARIABLE  DEFINITIONS 

C 

C  KS  CURRENT  AR  ORDER 

C  KT  CURRENT  MA  ORDER 

C  RX  AUTOCORRELATION  DATA  OF  INPUT  VN 

C  RY  AUTOCORRELATION  DATA  OF  OUTPUT  Y 

C  RXY  -  CROSSCORRELATION  DATA  OF  INPUT  AND  OUTPUT 

C  RXY  -  CROSSCORRELATION  DATA  OF  OUTPUT  AND  INPUT 

C  TXA  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 

C  OF  INPUT. 

C  TXB  VECTCR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 

C  OF  INPUT 

C  TYA  -  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 
C  OF  OUTPUT 

C  TYB  -  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 
C  OF  OUTPUT 

C  GA  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

C  OF  INPUT 

C  GB  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

C  OF  INPUT 

C  ZTA  -  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 
C  OF  OUTPUT 

C  ZTB  -  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 
C  OF  OUTPUT 

C  C  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXA 

CD-  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYA 
C  E  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GA 

C  F  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTA 

CP-  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXB 
C  Q  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYB 

C  R  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GB 

C  S  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTB 

C  A  -  VECTOR  UPDATED  AR  COEFFICIENTS  OF  ARMA  MODEL 

C  B  -  VECTOR  CONTAINING  UPDATED  MA  COEFFICIENTS  OF  ARMA  MODEL. 

C  TAU1  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 

C  ERROR  COEFFICIENTS. 

C  TAU2  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 
C  ERROR  COEFFICIENTS. 

C  TAU3  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 
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n  o  o 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ERROR  COEFFICIENTS. 

TAU4  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 
ERROR  COEFFICIENTS. 

XMU1  -  COEFFICIENT  OF  AR-TYPE  RECURSIVE  FORMULA 
YMU1  -  COEFFICIENT  OF  AR-TYPE  RECURSIVE  FORMULA 
MU2  -  COEFFICIENT  OF  AR-TYPE  RECURSIVE  FORMULA 
XMU3  -  COEFFICIENT  OF  AR-TYPE  RECURSIVE  FORMULA 
YMU3  -  COEFFICIENT  OF  AR-TYPE  RECURSIVE  FORMULA 
MU4  -  COEFFICIENT  OF  AR-TYPE  RECURSIVE  FORMULA 
MU5  -  COEFFICIENT  OF  AR-TYPE  RECURSIVE  FORMULA 
DET  -  DETERMINANT  OF  PREDICTION  ERROR  MATRIX  COMPOSED  OF 
VX.VY.VXY. 

ERR  -  ERROR  BETWEEN  REFERENCE  MODEL  OUTPUT  AND  LATTICE 
REALIZATION  OUTPUT. 

DIMENSION  RX(0:  2000) ,RY(0:  2000) ,RXY(0:  2000) ,RYX(0:  2000) ,A(0:  21) 
DIMENSION  B(0:  21) ,TXA( 0:  21) ,TYA( 0:  21),TXB(0:  21),TYB(0:  21),GA(0:  21) 
DIMENSION  GB(0: 21) ,ZTA(0: 21) ,ZTB(0: 21) ,C(0: 21) ,D(0: 21) ,E(0: 21) 
DIMENSION  P(0: 21) ,Q(0: 21),R(0: 21),S(0: 21),F(0:  21) 

REAL  MU5 ,MU2 ,MU4 


COMPUTE  VALUES  FOR  TAU1  THROUGH  TAU4 


T1S=0 

T2S=0 

T3S=0 

T4S=0 

KI=KS-1 

DO  10  1=0, KI 

T1S=T1S-RYX( I+1)*ZTA(KI-I) 
T2S=T2S+RY( I+1)*ZTA(KI-I) 
T3S=T3S+RY( I+1)*GA( I) 
T4S=T4S+RY( I+1)*ZTA( I) 

10  CONTINUE 

T1T=0 
T2T=0 
T3T=0 
T4T=0 

DO  20  J=0 ,KT 

T1T=T1T+RX( J+ 1 ) *ZTB ( KT- J) 
T2T=T2T-RXY( J+1)*ZTB( KT- J) 
T3T=T3T-RYX(KI-KT+1+J)*GB(J) 
T4T=T4T-RYX( KI -KT+ 1+J)*ZTB ( J) 
20  CONTINUE 

TAU1=T1S+T1T 

TAU2=T2S+T2T 

TAU3=T3S+T3T 

TAU4=T4S+T4T 

C 
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C  COMPUTE  VALUES  FOR  THE  REFLECTION  COEFFICIENTS 
C 

XMU1=-TAU1/VZ 

YMU1=-TAU2/VZ 

MU 2  =-VGZ/VZ 

DET  =VX*VY-VXY*VXY 

XMU3=-( VY*TAU 1 - VXY*TAU2 ) / DET 

YMU3=-( -VXY*TAU1+VX*TAU2)/DET 

MU4  =(  VGZ*TAU4-VZ*TAU3 )  /  ( VGZ*VGZ  -VG*VZ) 

MU 5  =MU4*MU2 
C 

C  COMPUTE  ARMA  COEFFICIENTS 

C 

DO  16  K=0 , KS 
C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 

16  CONTINUE 

DO  45  J=1 ,KS 
TXA( J)=C( J)+XMU1*F( KS - J) 

TYA( J)=D( J)+YMU1*F( KS - J) 

GA( J)=E( J- 1 )+MU2*F( J- 1 ) 

ZTA(J)=F(J)+XMU3*C(KS-J)+YMU3*D(KS-J)+MU4*E(J-1)+MU5*F(J-1) 
45  CONTINUE 

DO  31  K=0 , KT 
P(K)=TXB(K) 

Q(K)=TYB(K) 

S(K)=ZTB(K) 

R(K)=GB(K) 

31  CONTINUE 

DO  55  J=1 ,KT 

TXB(J)=P(J)+XMU1*S(KT+1-J) 

TYB(J)=Q(J)+YMU1*S(KT+1-J) 

GB( J)  =R( J)+MU2*S( J) 

ZTB (J)=S(J+1) +XMU 3*P( KT - J ) +YMU  3*Q( KT- J ) +MU4*R( J ) +MU5*S ( J ) 
55  CONTINUE 

C  WRITE  (*,176)  KS 

C  WRITE  (9,176)  KS 

176  FORMAT( 12) 

C  WRITE  (*,175)  ( ZTB(K) ,  K=0,KT) 

C  WRITE  (9,175)  (ZTB(K) ,  K=0,KT) 

175  FORMAT  (4( 1X.F10. 5) ) 

C 

C  UPDATE  ERRORS 

C 

650  FORMAT  (/'  S  UPDATE  ERROR  IS:  ' ,F15. 10) 

VX  =VX+XMU1*TAU1 
VY  =VY+YMU 1*TAU2 
VXY=VXY+XMU 1*TAU2 
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o  n  o 


VG  =VG+MU2*VGZ 

VZ=VZ+( XMU3*TAU1+YMU3*TAU2 )+MU4*TAU3+MU5*TAU4 
VGZ=TAU3+MU2*TAU4 
ERR=VY - ( VXY**2 ) / VX 
WRITE  (*,650)  ERR 
WRITE  (9,650)  ERR 
WRITE  (*,888)  VX.VY.VG.VZ 
888  FORMAT  (4(1X,F10.  6)) 

COMPUTE  MODEL  COEFFICIENTS 

DO  65  J=1 ,KS 

A( J)=TYA( J) -TXA( J)*VXY/VX 
65  CONTINUE 

DO  70  J=1 ,KT 

B(J)=TYB(J)-TXB(J)*VXY/VX 
70  CONTINUE 

B(0)=-VXY/VX 

RETURN 

END 
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APPENDIX  C.  SUBROUTINE  FOR  MAIN  PROGRAM 


SUBROUTINE  NEWMACKS.KT.VX.VY.VXY.VG.VZ.VGZ.TXA.TXB.TYA.TYB.GA.GB.Z 
&TA , ZTB , RX , RY , RXY , RYX , A , B ) 

C 

C  THIS  SUBROUTINE  COMPUTES  MA  PARAMETER  VALUES  FOR  AN  ARMA  MODEL  AS 
C  THE  MA  ORDER  INCREASES  BY  ONE. 

C 

C  VARIABLE  DEFINITIONS 

C 

C  KS  -  CURRENT  AR  ORDER 

C  KT  -  CURRENT  MA  ORDER 

C  RX  -  AUTOCORRELATION  DATA  OF  INPUT  VN 

C  RY  AUTOCORRELATION  DATA  OF  OUTPUT  Y 

C  RXY  CROSSCORRELATION  DATA  OF  INPUT  AND  OUTPUT 

C  RXY  -  CROSSCORRELATION  DATA  OF  OUTPUT  AND  INPUT 

C  TXA  -  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 

C  OF  INPUT. 

C  TXB  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 

C  OF  INPUT 

C  TYA  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 

C  OF  OUTPUT 

C  TYB  -  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 

C  OF  OUTPUT 

C  GA  -  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

C  OF  INPUT 

C  GB  -  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

C  OF  INPUT 

C  ZTA  -  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 
C  OF  OUTPUT 

C  ZTB  -  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

C  OF  OUTPUT 

C  C  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXA 

D  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYA 

C  E  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GA 

C  F  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTA 

CP-  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXB 
C  Q  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYB 

C  R  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GB 

C  S  -  ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTB 

C  A  VECTOR  UPDATED  AR  COEFFICIENTS  OF  ARMA  MODEL 

C  B  VECTOR  CONTAINING  UPDATED  MA  COEFFICIENTS  OF  ARMA  MODEL. 

C  TAU1P  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 
C  ERROR  COEFFICIENTS. 

C  TAU2P  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 
C  ERROR  COEFFICIENTS. 

C  TAU3P  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 
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C  ERROR  COEFFICIENTS. 

C  TAU4P  -  CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 
C  ERROR  COEFFICIENTS. 

C  XETA1  -  COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

C  YETA1  -  COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

C  ETA 2  -  COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

C  XETA3  -  COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

C  YETA3  -  COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

C  ETA4  -  COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

C  ETAS  -  COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

C  DET  -  DETERMINANT  OF  PREDICTION  ERROR  MATRIX  COMPOSED  OF 
C  VXjVY,VXY. 

C  ERR  -  ERROR  BETWEEN  REFERENCE  MODEL  OUTPUT  AND  LATTICE 
C 

DIMENSION  RX(0: 2000) ,RY(0: 2000) ,RXY(0: 2000),RYX(0:  2000) ,A(0:  21) 
DIMENSION  S( 0: 21),TXA(0: 21) ,TXB(0:  21),TYA(0:  21) ,TYB(0:  21),GA(0:  21) 
DIMENSION  GB ( 0 : 21) ,ZTA(0: 21) ,ZTB(0: 21) ,C( Or  21) ,D(0: 21) 

DIMENSION  E( 0: 21) ,F(0: 21) ,P(0: 21) ,Q(0: 21) ,R(0: 21),S(0: 21) 

C 

C  COMPUTE  VALUES  FOR  TAU1  PRIME  THROUGH  TAU4  PRIME 
C 

T1TP=0 
T2TP=0 
T3TP=0 
T4TP=0 
KJ=KT- 1 
DO  10  1=0, KJ 

T1TP=T1TP+RX( I+1)*GB(KJ-I) 

T2TP=T2TP-RXY( I+1)*GB(KJ-I) 

T3TP=T3TP+RX( I+1)*ZTB( I) 

T4TP=T4TP+RX( 1+1 )*GB( I ) 

10  CONTINUE 
T1SP=0 
T2SP=0 
T3SP=0 
T4SP=0 

DO  20  J=0 ,KS 

T1SP=T1SP-RYX( J+1)*GA(KS- J) 

T2SP=T2SP+RY(J+1)*GA(KS-J) 

T3SP=T3SP-RXY(KJ-KS+1+J)*ZTA( J) 

T4SP=T4SP-RXY(KJ-KS+1+J)*GA(J) 

20  CONTINUE 

TAU1P=T1TP+T1SP 

TAU2P=T2TP+T2SP 

TAU3P=T3TP+T3SP 

TAU4P=T4TP+T4SP 

C 

C  COMPUTE  VALUES  FOR  REFLECTION  COEFFICIENTS 

C 

XETA1=-TAU1P/VG 
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o  o  o 


YETA1=-TAU2P/VG 
ETA2  =-VGZ/VG 
DET  =VX* VY - VXY*VXY 
XETA3= - ( VY*TAU IP - VXY*TAU2P ) /DET 
YETA3=-(-VXY*TAUlP+VX*TAU2P)/DET 
ETA4  =( VGZ*TAU4P -VG*TAU3P ) / ( VGZ*VGZ - VG*VZ ) 

ETAS  =ETA4*ETA2 

COMPUTE  ARMA  PARAMETERS 

DO  16  K=0 ,KT 
P(K)=TXB(K) 

Q(K)=TYB(K) 

R(K)=GB(K) 

S(K)=ZTB(K) 

16  CONTINUE 

165  FORMAT  (4( 1X,F10. 5) ) 

DO  45  J=1 ,KT 

TXB(J)=P(J)+XETA1*R(KT-J) 

TYB( J)=Q( J)+YETA1*R(KT-J) 

GB(J)=R(J) +XETA  3*P ( KT - J ) +YETA3*Q( KT - J ) +ETA4*S ( J - 1 ) +ETA5*R ( J - 1 ) 
ZTB(J)=S(J-1)+ETA2*R(J-1) 

45  CONTINUE 

DO  41  K=0 , KS 
C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 

41  CONTINUE 

175  FORMAT  ( 5( 1X,F10. 5) ) 

DO  55  J=1 ,KS 

TXA( J)=C( J)+XETA1*E(KS+1-J) 

TYA( J)=D( J)+YETA1*E(KS+1-J) 

GA(J)=E(J+1)+XETA3*C(KS-J)+YETA3*D(KS-J)+ETA4*F(J)+ETA5*E(J) 
ZTA( J)=F( J)+ETA2*E( J) 

55  CONTINUE 

C 

C  UPDATE  ERRORS 

C 

VX=VX+XETA 1*TAU IP 
VY=VY+YETA1*TAU2P 
VXY =VXY +XETA 1  ,VT  AU  2  P 

VG=VG+(TAU1P*XETA3+TAU2P*YETA3)+ETA4*TAU3P+ETA5*TAU4P 

VZ=VZ+ETA2*VGZ 

VGZ=TAU3P+ETA2*TAU4P 

ERR=VY-(VXY**2)/VX 

WRITE  (*,66)  ERR 

WRITE  (9,66)  ERR 

66  FORMAT  (/'  T  UPDATE  ERROR  IS:  ’ ,F15.  10) 
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o  n  o 


WRITE  (*,889)  VX,VY,VG,VZ 
889  FORMAT  (4( 1X.F10. 6) ) 

COMPUTE  MODEL  COEFFICIENTS 

DO  65  J=1 ,KT 

B(J)=TYB(J)-TXB(J)*VXY/VX 
65  CONTINUE 

DO  70  J=1,KS 

A(J)=TYA(J)-TXA(J)*VXY/VX 
70  CONTINUE 

B(0)=-VXY/VX 

RETURN 

END 
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APPENDIX  D.  ADAPTIVE  LATTICE  ALGORITHM  PROGRAM 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  PROGRAM  CALCULATES  VALUES  OF  THE  LATTICE  COEFFI CENTS AND 
OUTPUT  OF  AN  ARMA  DIGITAL  LATTICE  FILTER  USING  AN  ADAPTIVE 
LATTICE  ALGORITHM. 

VARIABLE  DEFINITIONS 

X  -  ARRAY  OF  INPUT  DATA,  COMPUTER  GENERATED  WHITE  GAUSSIAN 
NOISE  WITH  UNIT  VARIANCE. 

AK  -  ARRAY  OF  LATTICE  COEFFICIENTS. 

R  -  ARRAY  OF  LATTICE  COEFFICIENTS. 

BO  -  TERMINAL  CONDITION  OF  LATTICE  REALIZATION. 

M  -  NUMBER  OF  LATTICE  STAGES  (EQUIVALENT  TO  ORDER  OF  ARMA 
MODEL). 

EXF  -  ARRAY  OF  FORWARD  PREDICTION  ERRORS  FOR  INPUT  X. 

EXB  -  ARRAY  OF  BACKWARD  PREDICTION  ERRORS  FOR  INPUT  X. 

EXBD  -  ARRAY  OF  DELAYED  BACKWARD  PREDICTION  ERRORS  FOR  INPUT  X. 

EYF  -  ARRAY  OF  FORWARD  PREDICTION  ERRORS  FOR  OUTPUT  Y. 

EYB  -  ARRAY  OF  BACKWARD  PREDICTION  ERRORS  FOR  OUTPUT  Y. 

EYBD  -  ARRAY  OF  DELAYED  BACKWARD  PREDICTION  ERRORS  FOR  OUTPUT  Y. 
ERROR  -  DIFFERENCE  BETWEEN  REFERENCE  MODEL  OUTPUT  AND  LATTICE 
REALIZATION  OUTPUT. 

YE  -  ARRAYS  CONTAINING  LATTICE  COEFFICIENT  VALUES  AT  EACH 
ITERATION. 

MU  -  CONVERGENCE  CONSTANT. 

RHO  *  WEIGHT  GIVEN  TO  CUURENT  POWER  LEVEL  AT  EACH  STAGE  OF  THE 
LATTICE  STRUCTURE. 

SIGK  -  POWER  LEVEL  USED  TO  NORMALIZE  CONVERGENCE  CONSTANT  WHEN 
UPDATING  AK  LATTICE  COEFFICIENTS. 

SIGR  -  POWER  LEVEL  USED  TO  NORMALIZE  CONVERGENCE  CONSTANT  WHEN 
UPDATING  R  LATTICE  COEFFICIENTS. 

SIGB  -  POWER  LEVEL  USED  TO  NORMALIZE  CONVERGENCE  CONSTANT  WHEN 
UPDATING  TERMINAL  CONDITION  BO. 

DIMENSION  EXF( 10) ,EXB( 10) ,EXBD( 10) ,EYF( 10) ,EYB( 10) ,EYBD( 10) ,R( 10) 
DIMENSION  AK(10),  X(9900),  ERR0R(9900),  V(12),  YE(6,9900) 

REAL  MU 
SIGK=1. 

SIGR=1. 

SIGB=1. 


INITIALIZE  ARRAYS 


DO  5  1=1,10 
EXF( I)=Q 
EXB( I)=0 
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VO  u  o  o  o  o  o 


EXBD(I)=0 
EYF( I)=0 
EYB(I)=0 
EYBD( I)=0 
R(I)=0 
AK( I)=0 
5  CONTINUE 

DO  6  1=1,9000 
X(I)=0 
ERROR(I)=0 
CONTINUE 

ENTER  VALUE  OF  THE  CONVERGENCE  CONSTANT  MU  AND  VALUE  OF  RHO. 

M=2 
N=300 

WRITE  (6,*)  'ENTER  MU' 

READ(6 ,*)  MU 
WRITE  (6,*)  ’ENTER  RHO’ 

READ  (6,*)  RHO 

GENERATE  WHITE  GAUSSIAN  NOISE 

ISIZE  =  N 
IX  =  152255 
I SORT  =  0 
MUL  =  2 

DO  7  K=  1, ISIZE 
CALL  SRND( IX, V, 12 ,MUL, ISQRT) 

XT=-6. 0 

DO  8  1=1,12 
8  XT=XT+V( I ) 

X(K)=XT 
7  CONTINUE 
C 

C  COMPUTE  OUTPUT  OF  REFERENCE  MODEL  FILTER  AND  LATTICE  STRUCTURE 
C  THEN  COMPUTE  THE  ERROR. 

C 

C  REFERENCE  MODEL 

Y3=0 

Y2=0 

Y1=0 

X3=0 

X2=0 

X1=0 

B0=1. 

DO  100  1=1, N 

C  YF=X( I)-0.  8*X1+1. 78*X2+0. 89*Yl-0. 25*Y2 
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YF=0. 5*X( I) -0. 4*X1+. 89*X2+0. 89*Yl-0.  25*Y2 
C  YF=X( I ) -2. 7*Xl+3. 21*X2-1. 59S*X3+1. 95*Y1-1.  62*Y2+0.  54*Y3 
C  YF=0.  5*X( I) -0. 95*X1+1. 33*X2-0. 979*X3+1. 69*Yl-0.  962*Y2+0.  2*Y3 
C  YF=X( I )+0.  2*Xl-0.  35*X2+1. 4*Yl-0.  85*Y2 
C  YF=0.  5*X( I) *0.  2*X1+Q.  445*X2+1.  0*Yl-0.  94*Y2 

C  Y3=Y2 

Y2=Y1 
Y1=YF 
C  X3=X2 

X2=X1 
X1=X(I) 

C 

C  LATTICE  FILTER 

C 

EXF( 1)=X( I) 

EXB( 1)=X( I) 

DO  10  K=1,M 

10  EXF(K+1)=EXF(K)+R(K)*EXBD(K) -AK(K)*EYBD(K) 

E YF( M+ 1 ) =B  0*EXF ( M+ 1 ) 

DO  20  K=1,M 

20  EYF(M+l-K)=EYF(M+2-K)+R(M+l-K)*EXBD(M+l-K) -AK(M+1-K)*EYBD(M+1-K) 

EYB( 1)=EYF( 1) 

DO  30  K=1,M-1 

EXB(K+1)=EXBD(K)+R(K)*EXF(K)-R(K)*EYF(K) 

30  EYB(K+1)=EYBD(K)+AK(K)*EYF(K) -AK(K)*£XF(K) 

YL=EYF( 1) 

ERROR( I)=YF-YL 
ERR=YF-YL 

CALL  UPDATE  ( R , AK , EYBD , EXBD , ERR , MU , M , S I GK , S I GR , B  0 , RHO ) 
CSB=EXF(M+1)*EXF(M+1) 

SIGB=RHO*SIGB+( l-RHO)*CSB 
B0=B0+( MU/SIGB )*ERR*EXF( M+l ) 

DO  40  K=1,M 
EXBD(K)=EXB(K) 

40  EYBD(K)=EYB(K) 

DO  50  J=1 ,M 
YE( J,I)=AK( J) 

50  YE( J+M, I)=R( J) 

202  FORMAT  (2( IX, F10. 6)) 

100  CONTINUE 

C 

C  PRINT  THE  ERROR  AND  VALUES  OF  THE  LATTICE  COEFFICIENTS. 

C 

WRITE  (*,200)  (ERROR(K) ,  K=1,N,10) 

WRITE  (9,200)  (ERROR(K) ,  K=1,N,10) 

200  FORMAT  (5( 1X.F10. 6)) 

WRITE  (9,209) 

209  FORMATC  ') 

WRITE  (*,201)  (R(K) ,  K=1,M) 
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WRITE  (*,201)  (AK(K) ,  K=1,M) 

WRITE  (9,201)  (R(K) ,  K=1,M) 

WRITE  (9,201)  (AK(K) ,  K=1,M) 

WRITE  (*,205)  BO 
WRITE  (9,205)  BO 
FORMAT  (F10.6) 

FORMAT  (5(1X,F10.  6)) 

CALL  PLOTTING  ROUTINES  TO  PLOT  ERROR  AND  LATTICE  COEFFICIENTS. 

CALL  PLOT  ( ERROR ,N) 

CALL  PLOT1  (YE,N) 

STOP 

END 

SUBROUTINE  WHICH  UPDATES  LATTICE  COEFFICIENTS. 

SUBROUTINE  UPDATE ( R, AK.EYBD , EXBD , ERR , MU , M , SIGK , SIGR , BO , RHO) 
DIMENSION  R( 10) ,AK( 10) ,EYBD( 10) ,EXBD( 10) 

REAL  MU 
CSK=0. 

CSR=0. 

DO  20  J=1,M 

C  SK=CSK+E YBD ( J ) *EYBD ( J ) *B  0**  2 
CSR=CSR+EXBD( J)*EXBD( J)*B0**2 
SIGK=RHO*SIGK+( l-RHO)*CSK 
S I GR=RHO*S I GR+( l-RHO)*CSR 
DO  10  J=1,M 

R(J)=R(J)+(MU/SIGR)*ERR*EXBD(J)*BO 

AK(J)=AK(J)-(MU/SIGK)*ERR*EYBD(J)*BO 

CONTINUE 

RETURN- 

END 

PLOTTING  ROUTINE  TO  PLOT  ERROR 

SUBROUTINE  PLOT(Y,N) 

DIMENSION  Y(N) ,X(9900) 

DO  10  J=1,N 
X(J)=J 
CALL  TEK618 
CALL  PRTPLT( 72 ,6) 

CALL  SHERPA( 'ADAPTIVE' , ’ A' ,3) 

CALL  RESET('ALL') 

CALL  PAGE( 8. 50,6.  0) 

CALL  HWROT('AUTO' ) 

CALL  XINTAX 

CALL  AREA2D(5.  0,3.  0) 

CALL  HEIGHT(0.  14) 

CALL  COMPLX 


on  o  ooni-*ooooon  non 


CALL  SHDCHRC 90. 0,1,0. 002,1) 

CALL  HEADINC ' LEARNING  CURVES? ' , 100 , 2. 0 , 1) 

CALL  XNAMEC ' ITERATIONS? ' , 100) 

CALL  YNAMEC ERROR?' ,100) 

C  CALLMESSAGC  ADAPTIVE  FILTER  ?' ,100,3. 0,-0. 8) 
CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF (0, 'SCALE' ,N,-3. 00, 'SCALE' ,3.  00) 

C  CALL  THKCRV(0.  02) 

CALL  CURVE(X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONE PL 
RETURN 
END 

SUBROUTINE  TO  PLOT  LATTICE  COEFFICIENTS 
SUBROUTINE  PLOT1  ( YE ,N) 

DIMENSION  YE(6,9900) ,X(9900) ,Y(9900) ,YD(9900) ,A(10) 

.  TRUE  VALUES  OF  THE  PARAMETERS 

A(l)=-0. 240719 
A(2)=0.  125 
A(3)=-0.  195719 
A(4)=0.  8900 
A(5)=0.  89 
DO  10  J=1,N 
0  X(J)=J 

CALL  TEK618 
CALL  PRTPLT(72,6) 

CALL  SHERPA( 'MENNECKE' , 'A' ,3) 

.  PRINT  SHERPA  FILE:  SHERPA  XXYYZZXX  SHGRAPH  A 

CALL  RESET('ALL') 

CALL  PAGE(8. 50,6. 0) 

CALL  HWROT( 'AUTO' ) 

CALL  XINTAX 
CALL  AREA2D(5.  0,3.  0) 

CALL  HEIGHT(0.  14) 

CALL  COMPLX 

CALL  SHDCHRC 90.  0,1,0.  002,1) 

CALL  HEADINC 'PARAMETERS?' ,100,2.  0,1) 

CALL  XNAMEC' ITERATIONS?' ,100) 

CALL  YNAMEC' MAGNITUDE?' ,100) 

CALL  MESSAGC 'ADAPTIVE  ARMA  LATTICE? ' , 100 ,3.  0 , -0.  8) 
CALL  THKFRMC0.03) 

CALL  FRAME 

CALL  GRAFCO,’ SCALE' ,N,-1.  00, 'SCALE' ,1.  0) 

CALL  THKCRVCO.  02) 

. TO  PLOT  ESTIMATES 

DO  20  K=1 ,4 
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30 


DO  30  J= 1 , N 
Y(J)=YE(K,J) 

CALL  CURVE(X,Y,N,0) 

20  CONTINUE 

C . TO  PLOT  TRUE  PARAMETERS 

DO  40  K=1 ,4 
DO  50  J=1,N 
50  Y(J)=A(K) 

CALL  DASH 

CALL  CURVE(X,Y,N,0) 

40  CONTINUE 

CALL  ENDPL(O) 

CALL  DONE PL 

RETURN 

END 


APPENDIX  E.  DERIVATION  OF  DIFFERENCE  EQUATION  FOR  TWO 

STAGE  ARMA  LATTICE  FILTER 

The  lattice  filter  is  described  by  the  expressions  for  the  forward  and  backward  pre¬ 
diction  errors,  equations  (3.6)  and  (3.7)  respectively, 

exA  (k)  =  x(k)  +  w2'  x(k  -  1)  -  w]  y{k  -  1) 
e*  (k)  =  exh  (k)  +  nj  e*  (k  -  1)  -  <  (k  -  1) 

4  (k)  =  x(k  -  1)  +  ivj  x(k)  -  wl  y{k)  (£  -  1 ) 

4,  (*)  =  3'(*  -  1)  ~  w,1  -t(A)  -f  iv]  y(k) 

eJ,  (*)  =  ^  (*)  +  »v4  %  (k-l)~  4t  (k  -  1) 

and  the  expression  for  the  filter  output, 

y(k)  =  (k)  +  wj  x(k  -  1)  -  wj  y(k  -  1)  (£-2) 

Substituting  for  e)x  {k)  in  equation  (E-2)  yields, 

y(k)  =  e*  (k)  +  nj  ex  (A  —  1)  —  nj  c*  (k  -  1)  +  wj  *(A  -  1)  -  wj  >(A  -  1)  (£  -  3) 

Substituting  for  c’}  {k  -  1)  and  (k  —  1)  in  equation  (E-3) 

y(k)  =  (k)  +  u'j  [  .v(A  -  2)  +  w’  ar(A  -  1)  -  w'  y(A  -  1)  ] 

-  M-v  [  r(A  -  2)  -  vv'  .r(A  -  1 )  +  w3  y(k  -  1 )  ]  {E  -  4) 

+  w4  x(k  -  1)  -  wj  y(k  -  1) 

Substituting  for  e}2  (A)  in  (E-4)  we  obtain 

>'(A)  =  eA  (A)  +  w\  e'l  (A  -  1)  -  wj2  ^  (A  -  1)  +  w42  x(k  -  2)  +  w4  w2  -t(A  -  I) 

— w4  w4  >'(A  —  1)  —  w2  y(A  —  2)  +  w2  w,1  x{k  -  1)  —  vv2  iv3!  >'(A  —  1)  (£  —  5) 

+w4  x(k  -  1 )  -  w3'  y(k  -  1 ) 

From  (E-l),  substitute  for  (A  -  1)  and  eix  (A  -  1)  in  (E-5)  to  obtain 
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(£-6) 


>•(*)  =  *JX  (k)  +  m-2  [  x(k  -  2)  +  wj  x(k  —  1)  —  wj  y(k  -  1)  ] 
-  w?  [  y(k  -  2)  -  wj  x{k  -  1)  +  wj  y{k  -  1)  ] 
4  wj  x(k  —  2)  4-  vv?  wj  x(k  —  1)  —  wj  wj  y(k  —  1) 
-ujyik  -  2)  +  wj  vv,1  x(k  -  I)  -  W32  W31  y(k  -  1) 
+  wj  x(k  —  1)  —  wj  y{k  -  1) 


Now,  from  (E-l),  substituting  for  e)  ( k )  in  (E-6)  yields 


>’(A)  =  jc(A)  +  wj  x{k  —  1)  —  w j  y(k  —  1)  +  wj  x(k  -  2)  +  wj  wj  x(k  —  1) 

-  vVj  wj  y(k  -  1)  -w  j  y{k  -  2)  +  wj  wj  x(k  -  1)  -  vv2  wj  y{k  -  1) 

4  vv2  x(k  -  2)  +  wj  wj  x(k  -  1)  — wj  wj  y(k  -  1)  -  wj  y(k  -  2) 

+  W3  w'i  *(*  -  1)  -  w3  w3  >'(k  ~  1)  +  wj  x(k  -  1)  -  wj  y(k  -  1) 


Grouping  terms  we  get, 

j-(A)  =  x(k)  +  (wj  +  wj  wj  +  wj  wj  +  wj  wj  +  wj  wj  +  wj)  x(k  —  1) 
4  (vv2  +  w42)  x(k  -  2) 

-  (vv,  +  W4  W2‘  +  W3  w."  +  w4  w4  +  w3  w3  +  w3 )  y(k  —  1) 

-  (wj  4-  wj) y(k  -  2) 


From  the  gradient  estimator  and  coefficient  update  equations  we  know  that  the  follow¬ 
ing  relationships  among  lattice  coefficients  are  true 


vv 


1  ! 


1  —  lv3 


2  2 
vv,  =  vv3 


W’-)  —  Wa 


{E- 9) 


Using  these  equalities  in  equation  (E-8),  we  obtain  the  final  expression  for  the  difference 
equation. 


y(k)  =  x(k)  +  2(vvj  4  wj  wj  4  wj  wj)  x{k  —  1)  +  2  wj  x(k  —  2) 
—  2 (wj  +  wj  wj  4-  wj  wj)y(A  —  1)  —  2  wj  y(k  —  2) 


(£-  10) 
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